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An experimental investigation of heat transfer effects 
on boundary layer separation in supersonic flow 


By G. E. GADD 


Aerodynamics Division, National Physical Laboratory, Teddington 


(Received 26 October 1956) 


SUMMARY 
Experiments have been done on the effects of heat transfer 
on wall-pressure distributions through separated regions with 
both laminar and turbulent boundary layers at a free-stream 
Mach number of about 3. ‘The temperature of the flat plate on 
which the boundary layer was formed could be varied from about 
—~35°C to +75°C. According to theory, this variation should 
have produced appreciable alterations at a laminar separation 
point in either the pressure or the pressure gradient, but no sign 
of this appeared in the overall pressure distributions, which, 
for laminar layers, remained unaffected by wall temperature. 
A possible explanation is given for this apparent discrepancy 
between theory and experiment. With turbulent layers, the 
variations in wall temperature did produce small changes in 
the pressure distributions. However, for most practical purposes 
such changes could be ignored. Hence the convenient conclusion 
is suggested that in supersonic separating flow with either a 
laminar or a turbulent boundary layer the pressure distributions 
are not significantly affected by moderate variations in wall 
temperature. 
1. INTRODUCTION 
Several theoretical papers (Cohen & Reshotko 1955a, 1955b; Gadd 
1952, 1956a; Illingworth 1954; Morduchow & Grape 1955; see .also 
Gadd 1956b) have been written on the problem of the effects on laminar 
separation of cooling or heating the wall. (By cooling is meant maintaining 
the temperature 7',, of the wall at less than the value 7',.. for which there 
is zero heat transfer between fluid and wall. 7°, is, of course, usually a 
little below the stagnation temperature.) All the theories agree in predicting 
that cooling the wall makes it more difficult to separate a laminar boundary 
layer, and that heating has the opposite effect. No theory has yet been 
developed for the effects of heat transfer on turbulent separation, but one 
would expect there to be effects qualitatively similar to those on laminar 
layers. However the effects ought to be of smaller magnitude, one would 
think, since, in the boundary layer upstream of the region of separation 
changes in wall temperature produce smaller proportional changes of 
velocity and density when the flow is turbulent than they do when it is 
laminar. With laminar layers, for an arbitrarily fixed pressure distribution, 
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it is typically found in the theories that the pressure coefficient* Cp, at 
the separation point roughly obeys the relation Cy, <(7,,/7,,.)-", where n 
is between 0-5 and 1, and 7’. is measured on the absolute scale. This 
considerable predicted effect of heat transfer is of great interest, because 
in many practical applications where boundary layer separation at supersonic 
speeds may occur the wall temperature will be much lower than the zero 
heat transfer value. Also, in certain wind tunnel investigations observations 


Shock- generating wedge 


















































~ 
Supersonic Alternative means of 
uinstream provoking separation 
— x5 
Dace no! for | id 
4 Spoiler Sf assage in plate for liquid 
i Bp A ee SES 
Duct for tunnel wall 
: ee boundary layer 
a L 
iner Xx 4 - aaa 
Sur porcs for plate arc sides’ LS Pipes fory] 
circulating 
liquid 
: x 
- x Se Sees Sek 
->Flow of li quic x 
! . sere 
~ _ a — et 








leading edge portion of non-conducting material 

















Plate as modified with insulated leading edge portion 
Figure 1. The flat plate on which the test boundary layer was formed. 


are made before the model and flow are in thermal equilibrium, so that 
errors may arise if the effects of heat transfer on separation are appreciable. 
Previous experimental work (for example Gadd, Holder & Regan 1954) 
on separation in supersonic flow has all been concerned with the insulated 
condition. Hence it was decided to investigate experimentally cases of 
laminar and turbulent separation with the wall heated or cooled. 

* The pressure coefficient C, 2(p—p.)/y-Mj p, where p = pressure, p, = free- 


stream pressure, y = ratio of specific heats, and ./, — free-stream Mach number. 
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‘The experimental arrangement was as shown in figures 1 and 2. The 
experiments were done at a Mach number of about 3. The test boundary 
layer was formed on a hollow flat plate through which hot or cold liquid 
could be circulated. Separation could be provoked either by a shock wave 
generated by a wedge held in the mainstream above the plate, or by a spoiler 
attached to the surface. With either arrangement separation occurs upstream 
of the agency provoking it. ‘he pressure distribution in the neighbourhood 
of separation is thus governed by the equilibrium between the thickening 
of the boundary layer and the associated deflection of the external flow 
from its original free-stream direction. Most of the theories for laminar 
layers mentioned above do not take account of this equilibrium process 
which will nearly always in practice govern the pressure distribution. 
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Figure 2. Pitot tube for boundary layer traverses. 
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Figure 3. Possible distributions of the pressure coefficient C,, for laminar separa- 
tions as a function of the Reynolds number R,, based on free-stream conditions 
and the distance x from the leading edge. 


However, one would expect the theoretical results for cases with arbitrary 
pressure distributions to have a qualitative relevance to the experimental 
laminar case. ‘Chey would certainly be relevant if it so happened that, for 
a given free-stream Mach number and a given Reynolds number at the 
upstream end of the compression region, the pressure distribution as far 
as the separation point always formed part of a single curve, independent 
of wall temperature, as shown in figure 3. Then, according to the theories, 
H 2 
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the pressure coefficient C,, at the separation point would increase as the 
wall temperature was reduced. It would seem natural in these circumstances 
for the pressure distribution downstream of separation to be of the shape 
shown in figure 3, and to vary with wall temperature as shown there, for the 
following reasons. When laminar separation occurs in practice, the pressure 
gradients downstream of separation fall off where the dead-air region becomes 
thick, and usually increase again further downstream because of transition 
to turbulent flow. Thus the pressure distribution shows a ‘laminar foot’ 
(Gadd, Holder & Regan 1954), with a fairly well defined ‘top’ at which 
the pressure coefficient C,,, may be readily determined. With zero heat 
transfer C,,, is roughly proportional to C,,, and at first sight there seems 
no reason why the ratio of Cp, to C,, should vary greatly with wall 
temperature. ‘Thus, for a fixed upstream pressure distribution as shown 
in figure 3, C),, as well as Cy,, would be expected to vary like 77", 
where n is between 0-5 and 1. However, more will be said on this point 
later. Meanwhile it must be borne in mind that it is far from evident that 
the pressure distribution upstream of the separation point ought in reality 
Pas Cooled wall 
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Figure 4. The variation with wall temperature of the pressure-coefficient distribution 
with laminar separation, according to Gadd (1956 a). 


to be independent of wall temperature, as shown in figure 3. According 
to Gadd (1956a) the pressure distribution varies with wall temperature 
not only downstream of the separation point, but upstream of it as well. 
Gadd’s theory is based on more realistic assumptions than the other theories 
mentioned inasmuch as it takes account of the equilibrium between the 
pressure gradients and the thickening of the boundary layer. It gives the 
result, which might appear to contradict the other theories for laminar 
layers, that the pressure coefficient at separation is unaffected by wall 
temperature, as shown in figure 4. However, a further result is that the 
pressure gradients at separation are sharper when the wall is cooled and 
more gradual when it is heated, since they vary as 7,,°*. If, with a given 
constant wall temperature, the pressure distribution could be changed 
from the form (1) of figure 4 to the form (3), without regard to the thickening 
of the boundary layer, separation would occur at a lower pressure with the 
form (3). Thus the theoretical result (Gadd 1956 a) that separation occurs 
at the same pressure with cooling despite sharper gradients is, in effect, 
consistent with the idea that cooling makes it more difficult to separate 
the boundary layer. 
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However, despite the qualitative agreement of all the theories, the 
experimental results for laminar layers appeared to contradict them. ‘The 
surprising result was obtained that the maximum pressure gradient at the 
upstream end of the laminar foot, the pressure coefficient C,,, at the top 
of the foot, and indeed the entire wall-pressure distribution, are all virtually 
unaffected by cooling or heating the wall over the range — 35° C to + 75°C. 
The wall temperature for zero heat transfer is always near 0°C, a little 
below the stagnation temperature (which is approximately atmospheric 
temperature). Hence the maximum cooling employed represents a pro- 


portional decrease of about § in absolute temperature, and the heating an 


increase of about }, so that the theories would certainly lead one to expect 
observable effects. Perhaps equally surprising is the result that there are 
bigger observable effects with turbulent boundary layers than with laminar 
ones. However, it may well be that there are real effects with laminar 
layers which are not apparent in the overall pressure distributions. Further 
discussion of the interpretation of the results is given in a later section. 


2. DESCRIPTION OF APPARATUS 

Figure 1 is a simplified diagram of the flat plate on which the 
test boundary layer was formed. ‘The plate completely spanned the 
2-6 in. x 1-5 in. blow-down tunnel (Gadd, Holder & Regan 1954) in which 
the investigations were made. Along the centre line of the plate there 
were pressure tappings and there were also thermocouples distributed over 
the surface as shown. ‘The thermocouple heads were soldered into 
= in.-diameter brass plugs which were forced into tightly fitting holes 
in the plate, and made flush with the surface. At the end of the investigations 
the plate was modified as in figure 1 (6) to have an insulated portion at the 
front. 

The liquid that was circulated through the + in.-deep passage in the 
plate was made to flow at about 5 to 10 ft./sec. Before entering the plate, 
it was pumped through a coil of copper pipe which served as a heat exchanger. 
The coil was placed in a large thermos flask so that it could be surrounded 
by liquid. When it was desired to cool the plate, methylated spirit was 
used for the circulating liquid and the liquid in the flask. Into the latter, 
lumps of solid carbon dioxide were dropped, and by this method the plate 
could be cooled to about — 35° C when low stagnation pressures were used 
for investigations on laminar layers near the front of the plate. For investi- 
gations on turbulent layers further along the plate, higher stagnation 
pressures were used to give the requisite higher Reynolds numbers. ‘This 
meant that the heat transfer was greater, and under these conditions the 
minimum temperature attainable was about —25°C. For heating the 
plate water was used as the circulating liquid and was heated by nearly- 
boiling water which was put in the flask. Plate temperatures of about 
75°C for laminar layers and 65° C for turbulent ones were attainable by 
this means. Heat conduction ensured that all metal parts of the plate 
surface (except perhaps those very close to the leading edge) took up a 
temperature which was roughly uniform. 
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Owing to the space required by the circulating liquid, the plate had to 
be thicker than the one used in previous investigations of interactions 
between shock waves and boundary layers (Gadd, Holder & Regan 1954). 
‘This increased thickness and the pipes ducting the liquid in and out of the 
plate made the blockage rather large, and at low stagnation pressures, like 
those used in the investigations of laminar flow, prevented the use of a 
shock-generating wedge of large angle. ‘Thus large regions of laminar 
separated flow were not readily obtainable by means of a wedge in the 
free stream, and the technique used in the former experiments of moving 
the wedge axially in order to move the pressure distribution past the pressure 
tappings was not entirely satisfactory. ‘This limitation was finally overcome 
by drilling extra pressure tappings. ‘The original spacing had been 3 in., 
and the later spacing was { in. in the region used for investigations of 
laminar flow and, finally, ; in. in the partly-insulated version of the plate 
shown in figure 1(6). ‘The more closely spaced holes made it possible to 
obtain satisfactory pressure distributions with laminar separations by 
using a fixed spoiler about 0-04 in. high glued to the surface at right angles 
to the flow and completely spanning the plate andtunnel. ‘This arrangement 
provoked separation about } in. to } in. upstream of the spoiler. 

A thin Pitot tube of external diameter 0-028 in. could be arranged to 
protrude through either of two of the static pressure tappings, one in the 
region where the laminar separations were made to oceur, and the other 
in the region of the turbulent separations. ‘The tube was bent over at the 
top and flattened at the tip to a height of 0-012 in. and a width of 0-034 in., 
as shown in figure 2. A micrometer underneath the working section moved 
it up and down so that pitot traverses of the boundary layer could be made. 


+ 


3. RESULTS AND DISCUSSION 
(a) Laminar separations 

Figures 5 to 9 show pressure distributions for laminar separations. 
‘The pressure p at the surface, made non-dimensional by dividing it by /p,, 
the constant value of p in the undisturbed region upstream of separation, 
is plotted against the distance x along the plate from the leading edge. 
The Reynolds numbers R, based on x and free-stream conditions are 
marked on the abscissa scales. It can be seen that there is no appreciable 
systematic effect of temperature in any of the figures, and mean curves 
have accordingly been drawn through the experimental points. ‘The curves 
all show the characteristic laminar foot, with moderate pressure gradients 
at the upstream end, a flatter top to the foot, and a steep increase of pressure 
at the downstream end. ‘This latter increase is due to transition to turbulent 
flow, which in the present experiments always occurred within the separated 
flow region when the flow was laminar at separation. 


Figures 5 to 7 were obtained with the all-metal plate as originally 
constructed, 1.e. they correspond to an approximately uniform wail tem- 
perature. With the plate modified as shown in figure 1(5) to have an 
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insulated portion near the leading edge, the wall temperature over the 
insulated portion always assumes the zero heat transfer value, near to 0° C. 
Hence when hot or cold liquid is circulated there is an abrupt step of wall 
temperature at the junction with the metal portion. Figures 8 and 9 
correspond to this condition. 

Figures 6 and 7 differ from the rest in certain respects. ‘The different 
figures correspond to distributions obtained on different occasions, the 
flat plate having been removed from the tunnel and put back again in 
between times. In figure 6 the laminar foot is rather short, the separated 
layer being more ready to turn turbulent than on other occasions. ‘This 
variation in transition behaviour may be due to a variation in the condition 
of the leading edge. An alternative possible cause is a variation in the 
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Figure 5. Pressure distributions with laminar separation induced by a wedge in the 
mainstream. Stagnation pressure 2 atmospheres. 


disturbance due to air leaking round the sides of the plate from underneath 
near the leading edge, although an attempt was always made to seal off 
this leakage. In figure 7 the non-dimensional pressure gradient 


re. ) 
‘ (6 


over the upstream part of the laminar foot is much greater than in the 
other figures. ‘This is probably due to disturbances emanating from the 
leading edge at the corners where it meets the side walls of the tunnel. 
Such disturbances would cross the plate centre line, on which the pressure 
tappings are situated, at about 2-1 in. from the leading edge at a Mach 
number of 3, and could well be responsible for distorting the pressure 
distributions a little upstream of this position. The other distributions 
are nearer the leading edge, and so are not affected. However, the distortions 
in figure 7, and the effects of the early transition in figure 6, should be 
independent of temperature, so that the two figures are admissible as 
evidence that temperature has little effect on the pressure distributions. 
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Figure 6. Pressure distributions with laminar separation induced by a wedge in the 
mainstream. Stagnation pressure 1:3 atmospheres. 
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Figure 7. Pressure distributions with laminar separation induced by a spoiler. 
Stagnation pressure 0-8 atmospheres. 
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Figure 8. Pressure distributions with laminar separation induced by a spoiler. 
J marks the junction between the insulated leading edge portion and the 
metal. Stagnation pressure 1:2 atmospheres. 
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Nevertheless, temperature does affect the boundary layer flow, as can 
be seen from figure 10. This shows Mach number profiles measured in 
the undisturbed boundary layer about 1-2 in. from the leading edge with 
the separation-provoking agency (wedge or spoiler) removed. ‘The Reynolds 
number R,, was 0-34 x 10° at the measurement position. ‘The Mach number 
was deduced from the static pressure at the wall and the readings of the 
pitot tube traversed through the boundary layer. In figure 10 the Mach 
numbers are plotted against the distances z, from the wall to the centre 
of the Pitot tube opening. Traverses made on different occasions did not 
repeat exactly, perhaps because of certain imperfections in the micrometer 
mechanism. For this reason the Mach number profiles for the different 
wall temperatures are first plotted separately showing the experimental 
14 
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Figure 9. Pressure distributions with laminar separation induced by a spoiler. 
J marks the junction between the insulated leading edge portion and the metal. 
Stagnation pressure 1:2 atmospheres. 


points, in figures 10(a), (6), and (c). ‘The different symbols for the points 
correspond to traverses made on different occasions. ‘The mean curves are 
superimposed in figure 10(d). Despite the scatter of the experimental 
points it is clear that there are real differences in the boundary layer profiles. 
The experimental points nearest the wall in figures 10 (a), (6), and (ce), lie 
above the mean curves, especially in the zero heat transfer and heated 
cases. ‘These points are obtained with the Pitot tube resting on the surface. 
The reading then corresponds to the Pitot pressure at a point further away 
from the wall than the centre of the tube opening. ‘Taylor (1938) showed 
that this displacement of the effective position of measurement decreases 


as the Reynolds number 
Ou 
a v(2 ) 
HB C3] wall 


increases, where d is the height of the Pitot at the tip, z is distance from the 
wall along the normal, and the other symbols have their usual significance. 
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For a given stagnation pressure and a given distance from the leading edge, 
the skin friction should be little affected by wall temperature (Howarth 
1953, pp. 418 to 425). The viscosity 1 is approximately proportional to 
absolute temperature, whilst the density p is inversely proportional to it, 
so that R, « 7,,* approximately, where 7’, is the absolute temperature at 
the wall. Hence, according to Taylor (1938), the effective displacement 
should be greatest with the heated wall, as is found experimentally. ‘The 
theoretical and experimental values for the magnitude of the displacement 


also agree reasonably. 
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Figure 10. Laminar Mach number profiles ; (a) with zero heat transfer (stagnation 

pressure 1-3 atmospheres, wall temperature about 6° C) ; (6) with a cooled 

wall (stagnation pressure 1-3 atmospheres, wall temperature about —34 C) ; 

(c) with a heated wall (stagnation pressure 1-3 atmospheres, wall temperature 

about 73°C); (d) the three curves superimposed (a cooled well, b zero heat 
transfer, c heated wall). 


The measured boundary layer thicknesses are about twice as great as 
the theoretical values for a flat plate boundary layer (Howarth 1953, p. 406). 
However, it appears from schlieren photographs that the boundary layer is 
probably not really as thick as the measurements indicate. The pitot tube 
used is rather large, with a tip height of about } the measured boundary 
layer thickness. Moreover, the tip is inclined downwards to the plate, as 
shown in figure 2, to ensure that it can make contact with the wall. Hence, 
since the tube is partially flattened at the tip so that its width is about three 
times its height, it must present a considerable disturbance to the flow, and 
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it is quite likely that it causes a local thickening of the boundary layer. 
Alternatively, the discrepancy between theory and the measured profiles 
may possibly be due to leading edge conditions. ‘The under surface of 
the plate near the leading edge (see figure 1) is inclined at about 18° to the 
upper surface, so that the pressures underneath must be much greater than 
those on the upper surface. Hence, since the leading edge shock must be 
slightly detached because of the inevitable slight bluntness of the edge and 
the presence of a boundary layer, it is possible that air leaks round from 
underneath at the leading edge, causing a small bubble of separation on the 
upper surface. ‘This would make the boundary layer thicker than the theory 
predicts (cf. Bradfield, de Coursin & Blumer 1954). 

The qualitative effects of heat transfer on the measured boundary layer 
profile shapes are similar to the theoretical ones, although the profile with 
zero heat transfer is closer to the heated wall profile than to that for the 
cooled wall, despite the fact that the temperature increase with heating is 
bigger than the drop with cooling. ‘This anomalous result is probably due 
to inaccuracies in the profile measurements. However, even if all the 
discrepancies between the theoretical and measured profiles are real, and 
due to some unknown experimental imperfections in the flow, it hardly 
seems likely that these imperfections can be such as to annul exactly all the 
theoretically predicted effects of heat transfer on separation. ‘Theory, as 
pointed out in $1, indicates that cooling the wall should in effect make the 
boundary layer more difficult to separate. According to Gadd (1956 a) the 
maximum pressure gradient at the upstream end of the laminar foot should 
be increased by cooling and reduced by heating in the ratio 7,,??. ‘Thus the 
pressure gradient for the cooled case should be about 1-8 times as great as 
that for the heated case with the temperatures of the present experiments. 
Such a large variation in gradient should be easily detectable experimentally, 
so that the present results definitely disprove Gadd’s predictions. However, 
his theory may in a sense be qualitatively correct inasmuch as, like the other 
theories, it predicts that cooling makes it more difficult, in effect, to separate 
the boundary layer. Indeed it seems inconceivable that all the theories 
could be wrong on this point. Hence it seems that the position of separation 
must occur at a higher pressure when the wall is cooled. This can still be 
the case even though the pressure distribution through the separated region 
remains substantially unaltered. ‘The precise position of the separation 
point cannot be determined merely by inspection of the overall pressure 
distribution, though it is known to be situated somewhere on the relatively 
steep upstream part of the laminar foot. It is possible, therefore, that the 
situation is as shown in figure 11, with separation occurring further upstream 
on the foot at the higher temperatures. In principle this could of course be 
verified experimentally using a surface pitot technique (Gadd, Holder & 
Regan 1954), but it is difficult to get accurate results by this method on the 
small scale of the present experiments. Hence it was not attempted with 
the present apparatus because considerable modifications would have been 
necessary and the blockage difficulties mentioned in §2 above would have 
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made the technique even more difficult than usual. It is hoped to make 
further investigations when a larger wind tunnel becomes available. 

At first sight it would seem natural for the pressures at the top of the 
laminar foot to be lower as in figure 3 if the pressure at separation was 
lower. If figure 11 represents the true state of affairs, the pressure gradients 
must fall off relatively more rapidly downstream of separation when the 
wall is cooled than they do when it is heated, and this requires explanation. 
A possible reason is as follows. 

When the wall is uniformly cooled the velocity profile of the upstream 
boundary layer has a ‘fuller’ shape than when the wall is insulated 
(cf. figure 10(d)). In the case of a plate as in figure 1 (4) with an insulated 
leading edge portion and a metal portion downstream which is cooled, the 
low-velocity part of the boundary layer rapidly assumes the fuller form on 
reaching the cooled part of the surface. This is because just downstream of 
the temperature discontinuity the viscosity at the wall is less than that in 
the middle of the boundary layer. Hence if the boundary layer profile 
remained the same as upstream there would be an unbalance of forces on 
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Figure 11. Suggested way in which heat transfer may in practice affect the position 
of the laminar separation point S. 


the fluid nearest the wall—the frictional force exerted by the wall tending 
to retard the inner layers of fluid would be outweighed by the force tending 
to drag them on exerted by the fluid further out. ‘This physical argument 
is confirmed by a mathematical solution which can be obtained for the flow 
near the wall near a temperature discontinuity. According to this solution 
the skin friction remains constant across the discontinuity if the viscosity 
is proportional to absolute temperature. Hence as the flow passes the 
discontinuity, the velocity gradient at the wall immediately assumes the 
value that it would have if the wall temperature over the insulated upstream 
portion were the same as that of the cooled part downstream. A little way 
from the wall immediately downstream of the discontinuity the velocity 
gradient remains approximately the same as if the wall were still insulated, 
so that the profile near the wall becomes curved, as shown in figure 12. 
However, the thickness of this region of curvature increases rapidly, varying 
as (v—.,)'%, where (x —.x,) is the distance downstream of the discontinuity. 
Hence the low-velocity part of the boundary layer profile very soon 
approaches the form that it would have if the wall were uniformly cooled. 
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This fuller upstream profile tends to make the boundary layer less easy 
to separate since the fluid fairly near the wall moves faster and is therefore 
less easily brought to rest. ‘This effect is reinforced by the fact that the 
density is increased by cooling, so that the momentum of the fluid is still 
further increased. On the other hand the viscosity near the wall is decreased 
by cooling. ‘his is a factor making for easier separation, because it is only 
viscosity that prevents the slowest moving fluid nearest the wall from being 
brought to rest immediately by the slightest pressure increase. However, 
the theories all agree that the combined effect of the fuller upstream profile 
and the greater density should outweigh the effect of the reduced viscosity, 
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Figure 12. Laminar velocity profiles near the wall near the juntion .J between the 
insulated and cooled portions of the wall. z distance normal to the wall, 


u velocity. 


and make the boundary layer more difficult to separate. Once separation 
has occurred, however, it seems quite possible that the viscosity effect may 
become dominant. The region downstream of separation is, as it were, 
further removed from the influence of the upstream profile. Hence it may 
quite well be that the pressure gradient which the boundary layer is able 
to withstand falls off more rapidly downstream of separation when the 
wall is cooled. ‘This could account for the fact that the overall pressure 
distribution remains unaltered, as in figure 11, despite the assumed increased 
pressure at separation. 

The apparatus sketched in figure 1(b), with an insulated leading edge 
portion, was intended to permit the effects of cooling or heating the wall 
in the separated region to be studied without making the upstream profiles 
fuller or less full. It was thought that with this arrangernent the pressure 
at separation, or at any rate that at the top of the laminar foot, would perhaps 
be reduced by cooling. However, as can be seen from figures 8 and 9, 
this was not found to be the case, probably because of the rapid adjustment 
in velocity profile that occurs at the temperature discontinuity, as described 


above. 
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One other point in connection with figures 5 to 9 deserves mention: 
the length of the laminar foot is virtually unaffected by wall temperature, 
and this means that the transition position is independent of heat transfer. 
It might perhaps have been expected that transition would be delayed by 
cooling. However, the 2-6 in. x 1:5 in.-tunnel in which the investigations 
were performed has a high turbulence level so that the stabilizing effect 
which cooling has on the laminar boundary layer oscillations, and the 
converse destabilizing effect of heating, may not much affect the position 
of transition. Alternatively it may be that since the length of the laminar 
foot is greater in terms of the upstream boundary layer thickness when the 
wall :s cooled (the thickness then being reduced), cooling does in fact have 
a stabilizing effect with regard to transition. 


(6) Turbulent separations 
Figure 13 shows pressure distributions for turbulent separations. The 
pressure p at the wall, divided by p,, the pressure in the undisturbed 
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Figure 13. Pressure distributions with turbulent separation induced by a wedge in 
the mainstream. Stagnation pressure 5:7 atmospheres. 


upstream boundary layer, is plotted against the distance x from the leading 
edge. ‘lhe Reynolds numbers R,, based on x and free-stream conditions 
are marked on the abscissa scales. Separation was provoked by a shock 
wave generated by a short 15 wedge in the mainstream. Since the wedge 
was short, the shock was followed closely by an expansion, so the peak 
pressures in figure 13 are not as high as those corresponding to the regular 
reflection of a shock of 15 flow deflection angle. ‘The pressure distributions 
have the characteristic shape for turbulent boundary layers with fairly 
large regions of separated flow (Gadd, Holder & Regan 1954). There is 
initial steep rise of pressure up to separation, downstream of which the 
pressure gradients fall off until the position where the shock strikes the 
boundary layer, and the pressure gradients then increase again until the 


peak pressure is reached. 
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Figure 14 shows Mach number traverses made at a distance of 3-8 in. 
from the leading edge in the undisturbed boundary layer, with the 
separation-provoking wedge removed. The Reynolds number R, was 
4-4 x 108 at the measurement position. The Mach number was deduced 
from the static pressure at the wall and the readings of the pitot tube shown 
in figure 2. In figure 14 the Mach numbers M are plotted against the 
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Figure 14. Turbulent Mach number profiles. Stagnation pressure 5-7 atmospheres. 
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Figure 15. The effect of Pitot size on measured Mach number profiles in a turbulent 
tunnel-wall boundary-layer. 


distances z, from the wall to the centre of the Pitot tube opening. The 
profiles are distorted near the wall due to pitot tube interference with the 
flow. This was established by traversing one of the thicker tunnel wall 
boundary layers with a small pitot tube and then with a large one whose 
dimensions and geometry relative to the tunnel wall boundary layer were 
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roughly similar to those of the tube used for the flat plate traverses of 
figure 14 relative to the flat plate boundary layer. ‘he tunnel wall results 
are shown in figure 15. ‘The results obtained with the large tube show too 
small a gradient near the wall, and too sharp a ‘shoulder’ a little further 
out. ‘There are evidently similar distortions in figure 14. Nevertheless 
the qualitative effects of heating or cooling are clear: heating the wall 
makes the boundary layer displacement thickness greater and reduces the 
velocities in the boundary layer. The slight reduction in free-stream Mach 
number in the heated cases is due to the static pressure at the wall (assumed 
to be constant across the boundary layer) being higher. ‘This is probably 
caused by the increased rate of growth of the boundary layer and the 
consequent increased deflection of the free-stream. 
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Figure 16. The definition of the kink pressure p; from the pressure distibution 
with a considerable extent of separated turbulent flow. 


The pressure distributions of figure 13 show that the upstream effect is 
greatest with the heated wall and smallest with the cooled wall. This is 
despite the fact that the maximum pressure at the downstream end of the 
distribution is somewhat reduced with the heated wall, and increased with 
the cooled wall. ‘The variation in maximum pressure is probably associated 
with the effects of the expansion emanating from the downstream end of 
the shock-generating wedge. When the wall is heated the boundary layer 
displacement thickness is greater, and the expansion can produce effects 
further.upstream, thus reducing the peak pressure. With a reduced peak 
pressure one would expect the shock to have a smaller upstream effect, 
other things being equal. However, the increased displacement thickness, 
besides enlarging the upstream region influenced by the expansion, also 
tends similarly to extend the influence of the shock further upstream, and 
this presumably more than counteracts the effect of the reduced peak 
pressure. 

The effects on the pressure coefficient at the separation point may be 
gauged by considering the kink pressure p,, defined as the intersection of 
the tangents of maximum and minimum slope, as in figure 16. It has been 
found in previous investigations (Gadd, Holder & Regan 1954) that p, is 
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close to the pressure at separation. From figure 13 it appears that the 
sink pressure ratio p,/P; 1S approximately equal to 2-2 for the cooled and 
nsulated cases, and to 2-1 for the heated case. ‘This variation may be merely 
due to experimental scatter, though by analogy with the theoretical pre- 


dictions for laminar layers it seems quite likely that there is some effect of 


heat transfer on the pressure coefficient for turbulent separation, such that 
neating reduces it and cooling increases it. However, much larger 


temperature differences than those used in the present experiments would 
I F I 


ve needed to produce very large effects. 


4. CONCLUSIONS 
he interesting result has been obtained that the pressure distribution 
it the wall with a laminar separation is approximately unaltered by wall 


temperature. It has to be remembered that the temperature range covered 


n the experiments was not very large, and in many practical cases there will 
’¢ a much more drastic degree of cooling. Furthermore, the experiments 
vere performed at one Mach number only (.V/ = 3), and with a high level 
of free-stream turbulence. ‘The latter factor is important since the position 
of transition largely determines the upstream effect for flows which are 


jJaminar at separation but which turn turbulent before reattachment. In 


the absence, however, of obvious reasons for doubting the general 
ipplicability of the present findings, it may be tentatively assumed that 
they are always approximately valid for cases of laminar separation in 
supersonic How with wall temperatures not too far from the zero heat 
transfer value, and with transition occurring 1n the separated layer. For 
such cases it would appear, therefore, that although as suggested by theory, 
there may well be a considerable effect of heat transfer on the pressure 
coefficient at separation, little sign of this appears in the overall pressur« 
listribution. For practical purposes the pressure distribution is usually 
f greater interest than the precise position of separation, so the experimental 
esult is very convenient. 
For turbulent layers, the pressure coefhcient at separation is probably 
fecreased by heating and increased by cooling, though the effect is not 
large. ‘he upstream effect produced by a given disturbance increases when 


the wall is heated and decreases when it is cooled. ‘This appears to be 


associated with the parallel effect of wall temperature on boundary layer 
displacement thickness, an effect which is only likely to be large when the 
vall temperature differs very greatly from the zero heat transfer value. 

Thus for estimating overall pressure distributions associated with 
youndary layer separation in supersonic flow when the wall is moderately 
1eated or cooled, it is probably sufficiently accurate in most cases to use 
jata obtained from experiments performed with zero heat transfer. 


The work described above was carried out in the Aerodynamics Division 


of the National Physical Laboratory, and is published by permission of the 


Jirector of the Laboratory. 


P.M. 
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The Saha equation and the adiabatic exponent in shock 
wave calculations 


"By RALPH A. ALPHER 
General Electric Research Laboratory, Schenectady, Nez York 


(Received 24 November 1956) 


lhe purpose of this note is to comment on the calculation of equilibriun 
vas dynamic parameters behind strong shock waves. ‘The writer has been 
ised by the appearance of a paper by Guman (1956) presenting 
neralized computing scheme for 1tonizing shock waves in monatomic 
gases. In that paper the reader is not cautioned about including excited 
states in the Saha equation for the computation of the degree of ionization 
chind shock fronts at appropriate temperatures and densities. ‘The same 
paper treats the adiabatic exponent y = c,,/c, as constant across strong 


i 


shocks when at the same time it is implied that the computing scheme is 
ot general validity. Hence, the unwary reader might attempt to apply 


the scheme in a regime where y is not only no longer constant but is no 
longer a useful quantity for characterizing the shock conditions. Othe 
thors (see, for example, Glass, Martin & Patterson (1953)) characteriz« 
s in which a shock has excited internal degrees of treedom in terms 
variable specific heat ratio when in fact one cannot use this quantity 

in calculating shock front conditions. 
he following comments on some of these points are presented with 
apologies to those readers who are familiar with so-called real gas effects; 


vever, the appearance of the material referred to above, as well as other 

mples in the literature, seems to require some comment. 

I.et us consider first the Saha equation. In a dilute electrically neutral 
gas in thermodynamic equilibrium at a temperature and density such that 
only neutral and singly ionized atoms are present, the concentrations 
(number density) of neutral atoms N,, of singly ionized atoms \,, and 
ot electrons .V,, are related by the Saha equation (see, for example, Aller 
(1%53), or Fowler (1955)) 

Ni N,/N, = 2(20m, RT/h?)??(B,(p, T)/BAp, T)jexp(—x,/RT), (1) 

vhere B, and B,, are the partition functions for atoms and ions respectively, 
factor 2 is the statistical weight for electrons, y, is the first ionization 
potential of the atom in question, and the other symbols have their usual 
meanings. At low temperatures and moderate densities, B,,/B, is 
adequately approximated by the ratio of the statistical weights of the 
ground states, g,,/g,. At either higher temperatures or lower densities, 
significant population of excited states occurs, and B,,/B, must be computed 
ore carefully. It may be noted that there appears to exist no unique and 
niversally accepted procedure for computing B, and B,,. Since the 


I2 
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partition functions for an isolated atom or ion diverge, one must invoke 
the presence of other neutral and charged particles in the gas to obtain a 
reasonable correction to the partition sums. ‘The recipes suggested by 
various investigators range widely in complexity (Bethe 1942; Fuchs, 
Lynch & Peierls 1942; Ecker & Weizel 1956; Meghreblian 1953; Unsold 
1948; Woolley 1955). Perhaps the simplest approach is to consider as 
bound to the atom or ion, and hence as contributing to bound-state 
partition sums, only those electrons whose classical orbits do not exceed 
half the mean separation of particles in the gas. Such an assumption, 
together with an adjustment in the statistical weight of the outermost 
orbit to allow for the fact that mean particle separation is not generally an 
integral number of Bohr radii, yields partition functions which agree 
satisfactorily with more detailed analyses which involve not only this 
‘pressure ionization’ but also electrostatic effects, modification of the 
more deeply bound energy levels, and correction of the partition function 
for free electrons (compare Fuchs, Lynch & Peierls (1942) and Bond 
(1954)). For example, the high temperature equation of state for argon 
computed by Bond (1954), using a * pressure ionization’ cut-off in the 
partition functions, agrees very well with shock wave measurements by 
Christian & Yarger (1955). We shall not attempt here to review the 
literature on the computation of partition functions, but rather would 
recommend to the reader who must deal with calculations of 1onization 


quilibrium the representative references which have been cited. 


Irgon Helium 

gm cm") T( WW) tO p(gm cm") T( BK) X Xp 
1-33 10 S$ OOO 1-040 10 50 OOO 1-11 
t +33: 5640 10 OOO 1-033 10 75 000 1-26 
2-66 10 12 000 1-029 5 10 ° 50 000 1°17 
2:66 - 10 16 OOO 1-029 5. 10 75 OOO 1-40 
4:00. 10 20 000 1-065 10> 50 000 1:20 
5-33.10 23 000 1-154 10 75 000 1-47 
§-33 x 10-3 25 000 1-225 10-4 25 000 1-01 

‘Table 1 


‘lo emphasize the effect of using only ground state statistical weights 
in the Saha equation, let us consider the calculation of the degree of 
ionization x, where 

a r Nf y tie if ’ >) 

x = N,,/(N,+Ny) = N./(N +N). (2) 
Let x, denote the ionization fraction calculated with corrected partition 
functions B,, and B,, and ~, the ionization computed using ground state 
statistical weights g,, and g,.. Then, if we write 


B, Bi, = (1 +f)g, 2 rIs (3) 
it follows that 
z,/%, = [1+(1—2)f} = (+f) (+23 Ay. (4) 
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based on partition functions computed by Bond (1954), and for helium, 
based on calculations by the author and colleagues. Both involve the 


‘This quantity is tabulated in table 1 at several values of p and 7 for argon, 


simple ‘pressure ionization’ correction. Clearly a deliberate decision is 
required before neglecting excited states in applying the Saha equation to 
strong shock calculations. 

Let us now turn to the use of the adiabatic exponent y in the 
Rankine-Hugoniot relations. ‘This index is usually introduced into the 
equations of conservation across a shock front by noting that the internal 
energy of an ideal gas is written 


I , 
E pd(1\p) = — P (>) 


5 Y l p” 
if p = kp’ and y isa constant. If (5) is valid across the shock, the Hugoniot 


condition may be written 


P» ’ Pa.) 
l(p,—D (— ps E.{ 1+ — E,( 14 Ee ) 
bs Pe Pi ~Py =) = Po E, ; Py E, 
(EF, — E,)y. (6) 
However, if the gas is vibrationally or electronically excited, dissociated, 
or otherwise changed in composition, or ionized by the shock, the internal 


enere\ 


ry is no longer a simple function of p and p as in (5), and the usual 


o(gm cm”) 10 10 10 10 10 
Helium, 25 000 K 1-46 £52 1-58 1-6 
1-65 1-65 1-66 1-66 
Helium, 30 000 K 1-22 146 1:17 1°26 
1-53 1-42 1-44 5() 
S000 K 17 1-18 1-2] | 3 
| -24 1-27 1-30 | () 
¢ I taken fro Cs yre (1955 ind Hirschtelde1 
( )48 
] 
S T1o t 7 ( mnt ) ry Fe t1o!1 I ak a ) 
ng NOSSID! (ne S S th strong Sn Gas 
Y. os 1¢ 
ntit Set B { x Velle1 (1940) \ 
ok p 7 
specific |! it rat s used by Gilmore 55) an 
Q a 55 ; 
sa we 1) | p Vk. (& 
tities require appropriate calculation of the high temperature 
( unic properties of the gas. Note that for an ideal gas y’ reduces 
isual Using (8) one may develop shock front relations involving 
1 } ] 1 , . ] 1 
for the unshocked gas (y, = y, ordinarily) and y, for the shocked gas 


’ which will formally resemble relations one might write dow1 


‘4 Yo} 
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with y, and y,. ‘Thus the usual Hugoniot relation can be written in terms 
of y’ in the form 
$[(P2— P:)(1/p1 + 1/p2)) = Bays - Ey. (9) 

Solution of this for p, p, will give an erroneous result if y, is taken as the 
specific heat ratio. ‘lo illustrate the errors possible, we present in table 2 
some comparisons of y and y’. 

Alternatively one may evaluate numerically the gas parameters behind 
a strong shock without introducing any special functions (see Resler, Lin 
& Kantrowitz (1952) and Romig (1956)). Evidently the principal use of 
the adiabatic index y in describing the thermodynamic state of a strongly 
shocked gas appears to be in calculating one of the sound velocities which 
mav be detined (see Bethe & ‘leller (1940)). 
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Buovant plumes in a moist atmosphere 


By B. R. MORTON 
Department of Mathematics, University College, London* 


(Received 24 September 1956) 


SUMMARY 
[his paper describes a simple model which can be used to 

investigate the transport of water vapour by thermal plumes 

in the atmosphere. For an approximate treatment of these 

plumes, it is assumed (as in a previous paper) that the vertical 

velocity, temperature and specific humidity are constant across 

the ascending column, and that the inflow velocity due to mixing 

at the edge of the plume is proportional to the vertical velocity 

within the plume. ‘lhe behaviour of the rising air is then 

investigated by means of equations representing the conservation 

of mass, momentum, heat and water vapour, and numerical 

solutions are obtained for representative cases. 

It is shown that in a stably stratified atmosphere the plume air 

will become saturated if the source is sufficiently strong, and that 

the height of this saturation level can be determined in terms of 

dimensionless parameters representative of the source and the 

environment. Above the saturation level the analysis is continued 

by taking into account both vapour and liquid phases of water, 

and an approximate treatment is given for the behaviour of small 

plume clouds. It is found that there is a critical size for these; 

smaller clouds remain characteristic of the plume, but larger ones 

extend upwards to heights which do not depend on the previous 

parameters, 

INTRODUCTION 

An account has been given recently by Morton, ‘Vaylor & ‘lurner (1956) 
of an approximate treatment for maintained buoyant plumes rising in a 
stratified environment under conditions such that the flow is turbulent. 
In this paper it was assumed that at each level the rate of entrainment at 
the edge of the plume is proportional to the mean vertical velocity at its 
centre. A similarity solution was then developed from equations repre- 
senting the conservation of mass, momentum and weight deficiency (relative 
to the ambient fluid at some specitied level, such as that of the plume source). 
The solution involves one undetermined constant, and for this a value 
was found from experiments in which a plume of dyed alcohol was released 
from the bottom of a tank containing stably stratified salt solution. Finally, 
the results of the theory were illustrated by a brief discussion of the ascent 
of thermal columns in the atmosphere. ‘This application to atmospheric 
convection will be extended in the present work. 


* Now at the Department of Mathematics, University of Manchester. 
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When a convective plume develops from a source in a moist atmosphere 
moisture will be entrained into the plume and carried up with the air. The 
resulting distribution of water vapour in the plume can be investigated by 
a simple extension of the methods of the earlier treatment, involving the 
introduction of a fourth conservation relation. Moreover, in the earth’s 
atmosphere the concentration of moisture is generally greatest near the 
ground (or sea) and decreases steadily with height, and thus at each leve! 
the plume air will be moister than its surroundings. Since the absolute 
temperature of the air in the plume also decreases with height under normal 
conditions, this air may become saturated at some level; a cloud will 
then appear in the upper part of the plume if suitable nuclei are present 
for condensation. 

his elementary model may possibly represent one mechanism for the 
formation of clouds, although certain reservations must be borne in mind 
if it is to be applied to the practical case. As the steady state rarely pertains 
to atmospheric phenomena such an approach can provide only a temporary 


ly varving condition. In addition, assuming that 


qaescription of a slow 
suitable plumes can form in the atmosphere, these rising columns will 
frequently be deflected by winds and under some conditions may lose 


their identity well below the heights predicted by the theory referred to 


above However, it may be noted that condensation does actually appeat 
to take pl n air which has risen trom lower levels. When atmospheric 
conditio sonabh I], clouds are sometimes seen above plumes 
rising from the chimneys of large industrial plants: a recent photograph 
of ( S¢ s] published by Scorer (1955), althoug! slightly 
a O Clouds Iso assoc ed t! 
n I ind’ thes vill scussed later lhus, 
th nsport « st vapour in buovant plumes discussed tn the present 
not SO pl I mong manv possible mechanisms for cloud 
toi led tl eather conditions are reasonably n and are 
nit 

Qe ( tv 
ot sou R ( 1 IS « ( univ | expected 
a Is plants and larg res d possibly 
ibe : irally o o sources is Mo! portant 
alti ses ut columns may be expected 
aby . 1 I ea sun-browned vegetatio 
W he s In each case the sourc ll 
SEN real the at transfer from some neighbourhood so that 
tl € | Mit se Wo may mu I Stronger than could be established 
by the source its In many cases these plumes will not penetrate to 
sufficiently great heights for condensation to be possible; they will then 


play an important but invisible part in the transfer of moisture in the 
atmosphere Another possible source of plumes in the atmosphere has 


been suggested by Batchelor (1954) and Priestley (1955), who have shown 


that convective plumes can grow spontaneously from small disturbances 
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in an unstably stratified layer. Although it seems likely that the effect of 
such a disturbance would be to start many plumes, and consequently t 
cause a general overturning of the laver, it is possible that heated air might 
continue to ascend above particularly vigorous regions in the layer and _ s« 
give rise to plumes which penetrate into the neutral or stably stratified air 
above. ‘These plumes will be at best quasi-steady because they have nx 
fixed points of origin and because they are likely to exhaust the layer locally ; 
but they may be of some importance in the transport of moisture upwards 
through the atmosphere. 

The analysis given previously (Morton, ‘Taylor & ‘Turner 1956, here- 
after referred to as paper 1) was developed for incompressible fluids, but 
an eXtension was made to cover convection 1n the atmosphere. ‘This was 
done by accepting the customary assumption that real convective motions 
in the earth’s atmosphere, described by velocities, potential temperatures 
and potential densities, are formally equivalent to motions in a similar 
but incompressible region described by velocities, ordinary temperatures 
and densities. While this substitution of potential temperature for 
temperature in the equations representing convective motion of an 
incompressible fluid compensates for the adiabatic nature of the pressure 
changes, the modified equations can be regarded as valid only within 


regions of limited vertical extent and for small total variations in potential 


temperature. For example, it can be shown that if the modified equations 
sed to investigate convection in the lowest two kilometres of the 
osphere the error in the solution is likely to be as large as 10 
spite of this, and although the ettects due to radiation have beet 
wholly neglected, it seems that the present type of analysis is likely t 
more accurate description of some aspects of atmospheric 
( ion than | Dec possib YEVIOUSI\ 


STEADY PLUMES IN A NEUTRALLY STABLE ENVIRONMENT 


ertain applications ti nvironment will be of uniform density 

tential density for gas), and this case will be considered first as 

( : xact solution to the approximate equations. ‘lhe motion above 

: source of heat tna still atmosphere at uniform potential temperatur« 
tery viindri polar coordinates (x,r), with the source at 

cin anc xis directec ly upwards. Let w be the vertical 
velocity, # the potential temperature and p the density inside the plume ; 
and let 4, and p, be the corresponding temperature and density 1n the 
ent air. lurther, take 0 4,(0) and p o,(O) as reference values fo1 


the system. Although the presence of water vapour in the atmospher 
Causes a small variation in the density and other.physical parameters, this 
etfect can be neglected, without introducing appreciable error, up to the 
height at which condensation begins. If it is now assumed that the velocity 
and buoyancy force are constant across the plume and zero outside it 
(i.e. a ‘top hat’ profile), the equations representing conservation of mass, 
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momentum and weight deficiency can be written (see paper I, equation 3) 


d 
Ty?) = 2xbu, (1) 
l 
Fy (bu) b?Be(0 — 4,), (2) 
Oe < 
y [b?u 32(4— 4y)] (), (3) 


) 


where 6 is the horizontal radius of the plume, 9 the coefficient of cubical 
g the acceleration due to gravity, and z the proportionality 


s 


‘xpansion, 
constant relating the inflow velocity at the edge of the plume to the vertical 
velocity within the plume. For boundary conditions it may be assumed 
that at the source the radius and the flux of momentum in the plume are 
zero, and that heat is released at a known and constant rate. Integration 
f equation (3) gives 


) 


o(4 — 6) F'/ 7, (+) 


h-ubd 


vhere F is the flux of buovaney. If p, is the density of the ambient 
tir at source level and c,, is the specific heat at constant pressure, the steady 
ite of output of heat from the source ts p,c,, fF Sg. The full solution of 


quations (1), (2) and (3) with the given boundary conditions can be written 


6% 5 /UxkF 


| 
6x \ Ts) 


Ji 


, SF /92F\ 
19(0 — H,) tain ( 


vd 67x\ 107 


Jt 


In paper I, x was determined tor a Gaussian profile since this titted closely 
vith the available experimental observations. ‘The equivalent value of z 
for the ‘top hat’ profile can be found by matching the vertical flux of mass 
ind momentum and the horizontal flow into the edge of the plume for the 
wo profiles, and is found to be x = 0-132. 
Phe distribution of moisture within the plume can now be investigated 
vv introducing a fourth equation representing the conservation of water 
vapour: 
; (b27qu) = 2bqy xu, (6) 
ax 
where q, the specitic humidity in the plume, is the mass of water vapour in 
init mass of moist air, and gy = g(x) refers to the surroundings. When 
relations (5) are used, equation (6) reduces to the form 


dq 5 5 a 
—s 4 
vith the solution 
5 : 2« { ) 
q = 3 constant + | x? 8q)(x) dx < (3) 
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One additional boundary condition is needed to determine the single 
disposable constant, and this is provided by the known rate of release of 
moisture from the source, which is proportional to E = [b?uq],,<». 

Although g,() may have any given functional form, it is sufficient for the 
present purposes to assume a linear profile for the specific humidity in the 
atmosphere and to take 





Jo(x) = A— Bx, (9) 

where 4 and B are profile constants. ‘Then equations (8) and (9) lead to 
; 5SE(9aF\-13 

a Ae Se 5/3, 0) 

J 1 — 2>Bx =a( to) \ (10) 


he first two terms of this expression depend on the moisture content of 
the atmosphere, while the third, which will generally be relatively small 
it greater heights, depends on the moisture released from the source. 
\loreover, equation (10) shows that all plumes from dry sources will have 
1 humidity distribution which is independent of the source strength. 

In order to find the level at which the plume becomes saturated it will 
be assumed, for the time being, that there is no dynamical etfect due to 
condensation; a subsequent correction can be applied above this height 

i if further information is required. ‘(he height at which saturation will 
ye attained can now be determined graphically on a tephigram by plotting 
‘elation (10) tor specific humidity together with the corresponding solution 
for the potential temperature of the plume air, 
5F (9aF\-13 
4=4,+ — ad lad (11) 
6z7aBe\ 107 


Vhis graphical determination of saturation heights is illustrated in 

gure 1 for a range of values of heat output (O = p,c,, F/ Bg = 0, 2°5 x 10°, 

+ 104, 1-5 x 10° cals per sec); each source is assumed to be dry (i.e. E = 0). 

Each of the heavily drawn curves (marked with its appropriate source 

strength) shows the variation in temperature of air rising in the buoyant 

plume up to the level at which saturation is reached, where it is terminated 

i it the heavy broken line which represents the dew point for plume air from 

iry sources of all strengths (see equation (10)). ‘The limiting case of zero 

urce strength represents the hypothetical ascent of a plume devoid of 

uoyancy, and it serves to detine a lower level below which cloud formation 

innot occur under the specified conditions. For purposes of reference, 

i second heavy broken line (to the left of the other curves) has been drawn 

to show the dew point of the environment, which is at uniform potential 
mperature. 

In figure 1 the scale of absolute temperature (7° C) is shown by 

ertical lines, that of potential temperature (@ C) by horizontal lines, 

that of pressure (p millibars) by continuous diagonal lines sloping to the 

ight and parallel with the ground level, and that of the saturation specific 

humidity (g¢,gm per kgm) by broken diagonal lines sloping to the left. 

It should be noted that the specific humidity used in the analysis has 
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measured the mass of water vapour in unit mass of moist air, in contrast 
with the usual meteorological unit of mass per kgm of moist air. For 
the purposes of the calculation the ambient air is assumed to be at uniforn } 
potential temperature, the temperature at the ground is taken as 20 C, 
and the profile of specitic humidity in the atmosphere is represented 
by relation (9) with the constants A = 11:-2x10-% gm per gm and 
B = 3-510 * gm per gm per cm. 
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T (°c) 





8 (°C) 
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Figure 1. ‘lhe ent of thermal plumes in a moist, neutrally stable atmosphere 
é natep m. ‘The continuous heavy curves show the tempera- 
f p e air above a ra »f source strengths (O); they are terminated 
level on a broken line showing the dew-point in all plumes 
t hows the variation with height lew-po 
vhic t ‘m potential temperature @ 20° ¢ 
{ \ af ed S | ) Curves that 1 
] | 1] ] 1 
ti I ti if plumes 1 re n vels where condens 1 : 
. a] ; a eee © 1 41 = cs 
I igh thes Vy hie a heights Devond the rang 
' i a 
pp mate theory l‘hus, in each case the formation of a cloud can ta 
+ lh] ] ¢ . o | . . 
plac suitable nuclei are present. However, the assumption of unifort 
tent ] ¢ —— . + @ . } t ~ = ] +; 
potential temperature for the atmosphere is rather too great a simplification, 
ind ti results shown hgure | are of value mainly because they illustrate 
the variation of specitic humidity with height lhe results of the analysis 
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may also be used to find the distribution throughout a buoyant plume of 
iny material that is being released from the source, or being entrained into 
the plume from the surrounding fluid. 


PLUMES IN A STABLY STRATIFIED ATMOSPHERE 


‘The atmosphere is commonly in a state of stable stratification: for this 
reason the treatment above will be extended in this section to the ascent 
of thermal plumes rising in a stably and uniformly stratified atmosphere. 
\lthough the method will be similar to that of the main analysis of paper | 
(pp. 6-10), the Gaussian profiles used there will be replaced now by * top hat’ 
protiles. This modification is necessary because two attempts to impose the 
normal profile in a treatment of buoyant plumes (Priestly & Ball 1955, and 
paper I) have led to errors in the solutions for the regions of negative 
buoyancy. Moreover, although the normal distribution profile probably 
provides a reasonable time average representation of the horizontal 
variation of vertical velocity and excess temperature relative to fixed axes, 
the mean profile relative to the plume axis is likely to be more nearly 
rectangular in form and of greater significance in finding the level at which 
the plume air saturates. 

It will be assumed, therefore, that the mean vertical velocity, excess 
temperature and excess specific humidity are constant across the plume 
and zero outside it, and can be represented by ‘top hat’ profiles for the 
approximate treatment of the present note. In deriving the equations 
it will be sufficient, for the present purposes, to make use of the previous 
ipproximations (see paper I, equations (7)) and, in addition, to ignore the 
effects of water vapour on the properties of air up to the level at which 
condensation begins. ‘The relations representing conservation of mass, 
momentum, heat, and water vapour can then be written 


1 
hd (b2u) = 2xbu, 
dx 
4 (b2u7) = Bgb?(6—4,), 
dx (12) 
d dé me 
— Bo 2 J — 3ah2 — 
= b*u( 6 — 4,)] oh u - 
sa i) — bey 2 40. 
dx q g, dx 


where g; = qo(0) provides a reference for qg. ‘To find the general nature 
f the convective motion, it is sufficient to restrict attention to the case 
in which both the potential temperature and the specific humidity in the 
indisturbed atmosphere vary linearly with height. ‘The environment is 
then characterized by the two parameters 


G = Pg d6,/dx and I= —(1/q,) dq/dx. 
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\s before, select the new variables 


; ; F 
| bu, i bu, F* = — = Bgb*u(6—6,), 





and VW a 
M* = — = #12 (13 
ae 41 
(in the present case F is the flux of buoyancy and M is related to the fluy 
of water vapour in the plume; p,c, Fy 8g and p,q, may be taken as ' 


the fluxes of heat and water vapour, respectively, from the source) 
Equations (12) can then be written in the simplified form, 
dit : dV} dF* dM* 

22V, 2F*W, = GW, — Wi.) =(14 
dx dx dx dx 


In thermal plumes of the type considered, there 1s zero flux of mass and 


momentum from the source, which 1s taken as a point, and a specified flus 

f heat and moisture. ‘Therefore, the boundary conditions at x 0 cat 
writtel 

| (), I] (), P =F VW Ve (15 

[Equations (14) and the conditions (15) are reduced to their simplest 


non-dimensional form by the transtormations 


2 Ay-12 f G a 
| ? | l FF “(| id 
I] 2 Ay V2 BAG 8g (1¢ 
[ ro hat 
he 7 ALF, Gm. 
It tollows that the most general form of the equations representing the 
scent of thermal plumes in a uniformly stratified atmosphere is 
dx dv ; df ; dm -_ _ 
a ee Oe ee et tp eee aos 
id the solutions to these equations must satisfy, at x, = 0, the boundary 


‘onditions 
GM, 


IF, 


vhere D is a dimensionless parameter. In spite of the changes that have 


dD, (18) 


(), w = 0), aoe ie m 


been made in equations (12) and relations (16) owing to the introductio1 
of the ‘top hat’ profile, the first three equations of (17) and first thre« 
conditions of (18) are precisely the same as the equivalent groups in paper I. 
Hence, the results of the numerical integration performed in paper I cai 


be carried over directly to the present case. Finally, a solution is required 


for the last equation of (17) satisfying the condition m = D at x, = 0. Th 
solution for f can now be used to infer that for m, and a set of values tor 7 
is shown in table 1 for the particular case of a dry source (i.e. D= 0). Te 
ybtain the equivalent solution for a source emitting water vapour, the 


relevant value of D should be added to each value in table 1. 
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All dry turbulent thermal! plumes are geometrically similar, although 
the scale of the convective motion depends on the parameters Ff, and G. 
In contrast, the effect of releasing moisture from a source in a damp 
atmosphere depends on Fy, G, J and M,; and hence D = GM)/1 Fy, may 
be regarded as a non-dimensional parameter which represents the dis- 
tribution of moisture in a thermal plume (or of the equivalent quantity in 
some corresponding system). A set of solutions for chosen values of the 


x Ty] 

0 0 1-0 00-1364 2:0 0-S487 
0-1 00-0003 1-1 0-1757 2:] 09632 
0:2 0-0019 1-2 00-2213 2-2 1-OS60 
03 0-0055 1-3 0:2736 2-3 1:-2173 
0-4 0-0199 1-4 0:3328 2:4 1-35608 
0-5 00-0215 1-5 00-3993 5 1-5045 
0-6 0:-0350 1-6 0:-4733 2-6 1 -6599 
0-7 0-0528 1-7 0-5550 Oe 1 -8227 
0-8 0-0753 1-8 0-6448 2:8 1-9912 
0-9 0-1031 1-9 00-7426 

lable 1. The numerical solution of the non-dimensional equation for the flux of 


a moist 


moisture in a thermal plume rising above a dry source of heat in 
atmosphere 


parameter D is shown graphically in figure 2; the non-dimensional quantity 


mw plotted is proportional to the value of (q¢—4q»)/q, on the axis of th 
plume. ‘The curve for D = 0 corresponds to a source emitting heat bu 
not moisture, and should apply to many natural sources; curves with 
dD 0 will apply, for example, to the plumes rising from certain industrial 
sources. ‘The following figures show orders of magnitude which might be 
ssociated with the column rising from the chimney of a large power station : 


zround temperature 24 C, constant decrease of temperature with height 


, 
nes 
gil 
| 
cure 2. The variation of excess specific humidity (proportional to m/z) with height 


(proportional to x,) for a thermal plume rising in a stably stratified atmosphere ; 
the various values of the dimensionless parameter D correspond to different 


rates of emission of water vapour from the source. 
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n the atmosphere 6-5 C per km, specific humidity at the ground 10 gi 
per kgm, gradient of specific humidity 3gm per kgm per km, heat output 
5 x 10° kW, water output (assuming half the heat wasted, and that 0-15 gm 
if water is produced in burning 1 gm of coal) 4-5 x 10* gm per sec. Under 
these conditions D = 3x 10°, although this value may be increased by 
water vapour drawn into the plume from the heat exchanger towers. 

To tind the approximate height at which the air in the plume first 
saturates it is necessary to consider particular cases. Figure 3, which is 
ilso plotted on a tephigram, shows a set of curves relating the temperature 
and height for the range of rates of heat output O = p,c,, Fy Bg = 41 « 10°, 
3-3 x 108 1 « 108, 2 10°, and 1:3 x 10!‘ kW from sources which emit no 
moisture (i.e. D = 0). Each plume is represented by a continuous curve 
corresponding with measurements of temperature and a broken curve 
showing the dew-point for air rising in the plume. In the lower parts ot 
) plume there is little effect due to stratification of the ambient air, so that 
these dew-point curves show individual behaviour depending on the source 
strength only in the uppermost part of each plume; hence, the lower part 
f each broken curve coincides with that drawn in figure 1 for the case of 
in unstratified environment, and these parts overlap in the diagram. ‘The 
dew-point curve for the environment has again been drawn (on the extrem« 
eft). Where the temperature and dew-point curves intersect, the plume 
tir becomes satura.cd with water vapour, although it can be seen that this 
happens only with the stronger sources (O 4-1 « 10° kW); when they 
lo not intersect, the two curves for a particular plume can be identified as 


they terminate at the same height. Conditions in the atmosphere have bee: 


taken as follows: ground temperature 20 C, ground specitic humidity 


li-2gm per kgm, G 1:09 x 10-4 c.g.s. units and / = 3:12x10-® c.g.s. 


inits (i.e. a decrease of 6-5 C and 3-5 gm per kgm in a kilometre of height). 
In determining these curves it has been assumed that air moved vertical 


\ 


indergoes dry adiabatic changes at all levels in the plume, so that no account 
has been taken of condensation and the consequent release of latent heat 
ibove the critical levels. It can be seen from figure 3 that for given 
itmospheric conditions the air in a thermal plume will ultimately saturate, 
provided that the source is sufficiently strong; these results give no 
nformation about conditions above the saturation level. When the source 
emits additional moisture the air will saturate at a lower level, although 
it can be inferred from figure 2 that this reduction in height will generally 
be small. 

The saturation levels predicted from these calculations depend on the 
ise of ‘ top hat’ profiles, and it can be shown that they are too high generally. 


If the excess temperature and excess specific humidity are actually greatest 
near the axis of the plume, the air will saturate first in this neighbourhood 
it some distance below the predicted height. Subsequently, as the plume 
tir rises, saturation will spread outwards from the axis; but it is difficult 
to estimate these effects adequately without a knowledge of the actual mean 
profiles of temperature and specific humidity in atmospheric plumes. It 
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nay also be remarked at this stage that some relatively poor approximations 
nave been used in the numerical examples; however, these are intended 
primarily as realistic illustrations of the model, and they are sufficient to 
show how to deal with particular cases of convection in the atmosphere. 


gure 
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3. A tephigram showing the behaviour of a family of plumes from sources of 
different strengths. For each source strength, a heavy continuous curve 


shows the plume temperature, and a broken curve the dew-point for plume 
air (these dew-point curves coincide except in the upper parts of each plume). 
Saturation occurs only when these two curves intersect below the height of 


the plume top. ‘The dew-point curve for the environment is again shown 


to the left. 


CLOUD FORMATION OVER CANE FIRES 


lhere is some difficulty in choosing a simple application for the 


calculi 


itions of the previous section, as the necessary observations are 


seldom all available. However, much of the important information has 


yeen 


obtained for one case, where cloud formation is known to be 


issociated with thermal plumes. It is apparently the practice to run a 
ire through the fields of sugar cane in Hawaii before harvesting, and 
zrowing cumulus clouds are often observed over the plumes of smoke 


rom 


W. A. 


these fires. The following information has been provided by 
Mordy of the Pineapple Research Institute at Hawaii. 


-M. 
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In a typical cane fire on a Hawatian sugar plantation, + x 104 square metres 
of cane is burnt in about halt an hour. On a field bearing 9 x 104 kgm of 
cane per acre, 1-8 x 10' kgm of leaves will burn in the fire. Half of this 
is water and will evaporate, and a further 4:5 x 10% kgm of water will be 
tormed by combustion, so that in all 1:35 10+kgm of water will be 
liberated per acre during the fire. Heat is also generated at the rate ot 
4:4 « 108 cals per kgm when moisture-free cane fibre burns, corresponding 
with 2-2 x 10® cals per kgm of the actual leaf material. Thus, during the 
half-hour of firing, water vapour will be released at the rate of about 
75 kgm per sec and heat at the rate of about 9-4 = 10° kW. In addition, 
tigure 4+ shows an atmospheric sounding taken in the early part of a morning 
on which cane-tire clouds were observed over Oahu. If no account is 
taken of conditions very close to the ground, the various parameters may 
be given the following approximate numerical values (in c.g.s. units) to 
represent the ascent of the thermal plume above such a burning cane field: 
Fi = Pee, G=POx«, f= 256x"*, My= 156x ond 
D = 0-025, 


re +. ‘The growth of plume clouds above sugar plantation fires at Oahu, Hawaii. 
An atmospheric sounding is shown with measurements of temperature (9) 
and dew-point (VY), and also a diagrammatic representation of the plume and 


cloud 


Before applying the analysis of the previous section, it will be necessary 

to make some allowance for the fact that a considerable area of the field 
ill be burning at any particular time. Although the determination of 
virtual sources has already been considered in paper I, there are several 
.dditional difficulties here since the moisture content of the air near the 
source must be known it a saturation level is to be predicted. ‘The effective 
size of the tire is unknown in this case, and it is difficult to estimate an 
appropriate area for the ground source because of the inflow of air near 
the ground and the emission of combustion gase ‘Table 2 shows a set 
of calculated results for the heights at which saturation occurs and to 
vhich the plume-top would rise if there were no condensation. The 


corresponding area of the source at ground level is also given for each of the 
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tual sources considered. ‘These results illustrate the general effect on 
the motion of the plume due to a large source; 1n the present case it seems 
hkely that saturation will occur at about 900 metres. However, no 
servation of cloud height was reported from Oahu, so that these 
calculations must be regarded as providing an illustration of the method 


ther than a test of the model. 





Ground area of 


depth of virtual Saturation height Height of plume-top 
s source 
irce (metres) (metres) (metres) 
(sq. metres) 
() () 1106 1482 
53 70 1064 1429 
106 280 1021 1376 
159 630 984 1323 
212 1120 947 1270 
265 1753 910 1217 
318 2530 873 1164 


! 











ple 2. The saturation level and height of the plume-top above the ground for the 
thermal plume over a burning cane field corresponding with several positions 
of the virtual source; the corresponding ground area of the source is shown 


in each case. 


\ similar calculation for a dry source (i.e. ) = 0) of small area indicates 
that the corresponding saturation level will be raised by no more than a 
few metres. According to this approximate theory, therefore, moisture 
released from the burning cane has little effect in producing saturation 
when this occurs at heights of about a kilometre. ‘Thus the formation of 
the cloud is due to the general conditions existing in the atmosphere and 

t to the emission of water from the source. 

An approximate profile for the plume has been marked in broken lines 

tigure + for the case of a small source at the ground, and a cloud has 
been sketched in to indicate the appropriate condensation level. As water 
condenses the temperature in the cloud will rise, and the consequent 
ncrease in the buoyancy forces will cause the uppermost parts of the cloud 
to extend above the original top of the plume. ‘This effect will be considered 


the tollowing section. 


"THE GROWTH OF CLOUDS IN THERMAL PLUMES 

In the previous sections an approximate theory of thermal plumes in the 
mosphere has been extended to show the height at which the rising ait 
become saturated. Above this level it may be assumed that a cloud 

\] develop if suitable condensation nuclei are present: Fora more complete 
nvestigation of the clouds which grow in thermal plumes, it will be necessary 
allow for the way in which condensation and evaporation effects vary 
roughout the condensed cloud, and this will require improvements in 
he model. However, some information can be obtained for the earlier 


K2 
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stages of the development of such plume clouds by a simple extension of 
the methods used already. 

For this purpose it will be assumed that there is a steady (or slowly 
varying) cloud in the upper regions of a thermal plume, and a moditied 
set of conservation relations will be applied above the condensation level. 
Although the cloud is unlikely to occupy the full width of the plume it 
will be convenient to take an averaged distribution of liquid water with a 
rectangular profile of the same width as the plume. It will be assumed 
also that the condensed water remains in very small droplets which move 
with the neighbouring air, and that only a small proportion of the water 
vapour condenses; consequently, both the relative motion of the drops 
and plume air and the small reduction in volume caused by the condensation 
will be neglected. 

‘The equations representing the conservation of mass, momentum, heat, 
and water (including both gaseous and liquid phases), can now be extended 
from their previous forms (12) to allow for the liberation of latent heat on 
condensation and the presence of liquid water. As the mass of liquid water 
per unit mass of the cloud mixture ts small, they may be written approximately 


is 


i 
(b2u) = 2bau, 

dx 

d 

oe (A277) h2 32(4 6.) heg y, 

ax 
; , j (iY) 
( Be ho - ' IL ks 6 es 
~ [b?uBg(8— 4)] b?u( 8g dé, dx) + —— = (b?uc), 
aX Cy dx 


d 


d 
(h2uq) 4 (b2uc) Dhxau n> 
dx 4 \ : 


d 
where po 1s the mass of liquid water droplets per unit volume, Z is the latent 
heat of vaporization for water, and c,, is the specific heat at constant pressure 
for air. ‘lhe variation of absolute temperature in small clouds will be small, 
so that L and ¢,, can be taken as constants for the present purposes. 

One more relation is required before a solution can be obtained, and this 
represents the condition for saturation within the plume; it can be written 
.n the form 
Gs = 9.(%, F), (20) 


where g, is the specific humidity required to saturate air at potential 
temperature @ and height x. It has already been pointed out that the 
variation of temperature within the cloud will be small provided that the 
treatment is restricted to the early stages of cloud development; under 
these conditions g, can be expressed by a simple approximate algebraic 
form. 

Just as equations (12) were reduced to the dimensionless set (17), the 


same transformations, with the additional relations 


We=K and K=a'*g"F,h, 
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can be used now to reduce equations (19) to dimensionless form: 
dw 
dx, 
dv! 
dx, z 
: (21) 
df dk 
dx, 
dm _ dk 
— +N 
dx, dx, 
BLc, and N = G_gq,1 are two additional non-dimensional 
The first depends only 








where = 
parameters needed to specify the cloud system. 
on the physical constants of air and water and will vary little from one 
plume cloud to another, while the second parameter provides a measure 
for the relative stratification of temperature (or of density) and moisture. 

In order to obtain a suitable form of equation (20) for use in numerical 
calculation, it will be assumed that the saturation water-vapour pressure 
is linearly related with the temperature over the range of a few degrees 
centigrade which may occur within the plume cloud. Although such a 
lationship is not strictly correct, it will be sufficient for the present 
made for the variation 


approximate treatment. When due allowance is 
of pressure with height, equation (20) can be written as a_ relation 


between q. and @: 


(G.— Yo) 4 V(x)(9— 6,) + X.(¥), 
hich the coetticients .Y, and _Y, are functions of the height x, and will 
pend on the particular case considered. Under the transformations used 
dy this reduces te 
m V(x) f+ VXo(x,)e, (22) 
A, and \, are appropriate functions. 
While a numerical solution can be found without difficulty for 
quations (21) and (22) in any particular case, a general solution cannot 
given because of the occurrence of two unknown parameters (” and .\) 
quations (21) and light variations of the coefficient functions in 
tions (22), and also because the boundary conditions at the cloud bas« 


ery sensitive to eht at which saturation occurs first. ‘hese 
v conditions require that there must be continuity in the flux of 
omentum, heat, and water vapour, across the saturation level, and 
of liquid water. When these are transformed into conditions on 
dk, it is apparent that the nature of the solutions will depend 


much on the height x, at which saturation o¢curs first in the plume 


f 
» Wi an 


see table 1, paper 1). 

here are two kinds of behaviour in the plume cloud, and to illustrate 
se it is necessary to describe the results of two calculations. In each 
se it will be assumed that there is a thermal plume rising above a source 


t small area and strength 10°kW in an atmosphere with ground temperature 
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20) C and a uniform decrease with height of 6°5° C per km. [n the nrst 
case the specific humidity will be taken as 12-23 gm per kgm at the grou: 
ind with a constant rate of decrease with height of 3-5 gm per kgm per 

It may be shown that saturation occurs at vy, = 2:5 (1.e. at a height of ab 
1320 metres); above this level, a numerical analysis must be perfor 

in order to obtain a new solution for equations (21) and (22). ‘This 1 

be described conveniently in terms of the following non-dimensional ratios 
wv, which is proportional to the plume radius b; 7? «w, which 1s proportic 

to the vertical velocity uw in the plume; f z, which is proportional to t 
excess of temperature (9—4,); mw, which is proportional to the excess 
specific humidity (q¢—q,)); and kw, which is proportional to the concen- 
tration o of liquid water in the plume cloud. “Vhe physical variables 


then be recovered from the relations 


v= 1-341 FG **x,, 
bh = 0:354FlIG 3 sy 7 
u = 1-896 F) 4G! 52? /a, 
39(0 — 0,) = 1-341 FLAG 5f /ze, 
(g—q)/q, = 1341 FG 3°Im zw, 
go = 1:341 F} 4G 8k/e. 

Figure 5 illustrates the case of a plume in which saturation occurs 
the level corresponding with x, = 2:5. ‘lhe continuous curves show the 
solutions for equations (17) below x, = 2:5, and for equations (21) and (22) 
above that level. In order that the effects of condensation may be shown 
clearly, the solutions to equations (17) (i.e. for the unsaturated plume) 
are continued above the saturation level by means of broken lines. Although 
the saturation level is some 160 metres below the top of the undisturbed 
plume, it can be seen that the only etfects due to the presence of the plume 
cloud are that the cloud top rises an extra forty metres, and that the internal 
temperature increases slowly with height (instead of decreasing) although 
the air inside the cloud still remains colder than its surroundings. 

The second case, illustrated in figure 6, is for a plume of the same 
strength ascending in an atmosphere which is similar to that of the previous 
case except that the ground specific humidity is now taken as 12°84 gm per 
kgm; the corresponding saturation level is x, = 2-1 (a height of about 
1110 metres). It can be seen that in this case the cloud will dominate 
the behaviour of the upper part of the plume. Within the plume cloud the 
excess temperature (and hence the buoyancy) increases with height, so 
that the cloud extends upwards until its uppermost parts reach an inversion 
strong enough to prevent further ascent. 

It is now possible to make certain general statements about plume 
clouds. If saturation occurs first near the top of the plume a small cloud 
can exist in the uppermost part of the plume without producing very much 
effect on the air flow. For lower saturation levels almost the only ettect 


is that the steady clouds will be correspondingly taller; in all these cases 
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Figure 5, The ascent of a thermal plume for which the saturation level corresponds 
tox, = 2:5. Below this level the curves illustrate the ascent of an unsaturated 
plume, while above the saturation level these solutions are continued as 
broken curves and the actual behaviour of the air in the cloud is shown by 
continuous curves. ‘The non-dimensional variables plotted are : w/v X plume 
radius; 7°/«w & vertical velocity; f/w X% excess temperature; m/z X% excess 
specific humidity; and k/w < liquid water concentration within the plume. 
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F:gure 6, The ascent of a thermal plume for which the saturation level corresponds 


to x, — 2:1. The variables plotted are as in figure 5 
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the clouds may be regarded as controlled by the plumes. However, when 
the saturation level occurs below a critical height, which appears to be a 
little above a, = 2:2 (i.e. 0-03 F)*G 3* metres), the behaviour of the flow 
in the upper part of the plume is changed completely. ‘The vertical extent 
of the plume cloud is no longer determined by the parameters characterizing 
the plume, but must depend on factors which have been omitted in the 
present simple model. ‘These will include variations in the structure of 
the environment, such as strong inversions, and possibly the more compli- 
cated pattern of flow which will occur in real clouds. It may be noted in 
passing that this prediction of a critical size for clouds may also be of 
importance when the clouds are formed in other ways, for example in 
lee Waves. 

‘The present treatment has been developed for a steady (mean) state, 
and it cannot be used directly to provide information about the ascent of 
plumes above sources of steadily increasing strength, and about the 
consequent process of cloud growth that may take place. However, the 
conclusions can be applied to give a qualitative idea of what may happen 
in slowly varying systems. In addition, while the results obtained should 
be reasonable for clouds which are controlled by their parent plumes, they 
are unlikely to be satisfactory for clouds which dominate the upper parts 
of their plumes. For a description of these larger clouds, it will almost 
certainly be necessary to make more detailed allowance tor the liberation 
and absorption of latent heat. It is clear that the liberation of heat in the 
lower and inner parts of this section of the plume, and the absorption in 
the upper and outer parts, will cause a considerable modification of the 
temperature distribution which has been assumed. ‘Uhere will be a tendency 
to increase the upward plume velocity at the bottom of the cloud and in 


the centre, and to reduce it at the t 


p and around the sides. When a cloud 


grows past the critical size it may be taken that the increased buoyancy 


forces in the centre and the decreased forces at the edge dominate the effects 
due to the plume. ‘lhe cloud is then no longer properly a part of the plume, 
and some asymmetry in its structure or some wind etkect may cause it to 
drift clear, or it may persist when the plume dies down. From the need 
for contin s apy nt that the central upwards motion must now be 
balanced by downdrafts at the sides, and though the motion in the cloud 
is turbulent it will possess to some degree the elements of a spherical vortex 
(not all of which will be marked by condensed water droplets). A completely 
new model will be necessary in order to investigate the behaviour of these 


deve lope d ( louds 
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SUMMARY 

In this paper, two-dimensional! problems of the diffraction and 
stability of shock waves are investigated using an approximate 
theory in which disturbances to the flow are treated as a wave 
propagation on the shocks. ‘Vhese waves carry changes in the 
slope and the Mach number of the shock. ‘lhe equations 
governing the wave propagation are analogous in every way to 
the non-linear equations for plane waves in gas dynamics, and their 
solutions can be deduced by the same mathematical techniques. 
of the waves is found to be an 


Since the propagation speed 
increasing function of Mach number, waves carrving an increase 
in Mach number will eventually break and form what we may 


call a ‘shock’ corresponding to the breaking of a compression 


wave into a shock in the ordinary plane wave case. Such a 


‘shock’ moving on the shock is called a shock shock. ‘Vhe shock 


shock is a discontinuity in Mach number and shock slope, and it 
must be fitted in to satisfy the appropriate relations between these 
} 


discontinuities and its speed. ‘he waves moving on the shocl 
are interpreted as the trace of cylindrical sound waves in the flow 
behind the shock. In particulat a shock-shock is the trace of 

enuine shock in the flow behind, and thus corresponds to Mach 


reflection. 

‘The general theory of the wave propagation 1s set out in § 2. 
Lhe subsequent sections contain ipplicat ons of the theory te 
specific problems, including the motion of a shock along a curv 
wall, ditfraction by a wedge, stability of plane shocks and th 
] shor K. 


instabilitv of a converging evlindric 

1. INTRODUCTION 
In this paper a relatively simple approximate method is developed for 
treating problems of the diffraction and stability of shock waves. ‘The 
theory can be formulated without reference to any specific problem and 


is convenient to give the basic ideas before discussing the applications. 
the 


1t 
it 


Only two-dimensional problems are considered in this first part; 


extension to other cases is given in Part II. 
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We start by considering the set of curves formed by the successive 
Sositions of a curved shock as it moves forward through a uniform medium, 
ind we introduce the orthogonal trajectories of this set of curves. ‘hes 
yrthogonal trajectories will be called ‘rays’. In figure 1 the positions of a 
shock moving trom left to right are shown as full lines and the rays are shown 
is broken lines. ‘This network of shock positions and rays may be used as 
i basis for orthogonal coordinates in the plane, and accordingly coordinates 
(z, 8) are introduced such that the shock positions are the curves x = constant, 

¢ 


ind the ravs are 8 = constant. A suitable z coordinate would be the time 1 


it which the shock occupies that position, but we modify this slightly and 





Figure 1. Sketch showing the successive positions of a curved shock ; the full and 
broken lines represent the shock positions and the rays, respectively. 


take « = ayt where a, is the sound speed in the uniform gas ahead of the 
shock. ‘Then, the distance along a ray between the shock positions given 
by x and ~+dz is M(x, 8)dx where M is the Mach number of the shock 
at (a, 8). If we let A(z, 8)dS be the corresponding distance between the 
rays 8 and 8+ d8, then it can be shown that, tor purely geometrical reasons, 
M and A must satisfy the ditferential relation 


é é- eA\ = = , 
sit ia) * 2 AE) sn ”) 


An elementary proof of this will be given in § 2 but froma more sophisticated 
point of view it is the condition for the space to be flat; the curvature 
tensor can be expressed in terms of Mand A, and (1) is the only component 
which is not identically zero. ‘his alternative derivation is given in some 
detail in Part II since it is the neatest way of obtaining the relations 
corresponding to (1) for three dimensions. 

Now, if we can tind a second relation between .W/ and 4A, we have an 
explicit equation for the Mach number of the shock as a function of (2, 8); 
from this function, the shock position can be determined for all times. 
The second relation must come trom the dynamics of the motion and 
strictly requires a solution of the equations of motion for the flow behind 
the shock, subject to the Rankine-Hugoniot relations across the shock and 
boundary conditions at solid walls, etc. Of course this is the original 
problem. But the above approach suggests a simple approximate procedure. 


© 
& 
5 








A new approach to problems of shock dynamics 147 


o some extent, the propagation of the shock between any two neighbouring 
avs can be treated as if the rays were solid walls. ‘Vhis would be exactly 
true if the rays were particle paths, but the most we can say is that immediately 
behind the shock the particles move normal to the shock, 1.e. in the rav 
lirection. However, we assume that the later divergence of the rays and 
the particle paths 1s not important and accept this similarity to propagation 

achannel. Now, for a shock moving down a channel, if the modifications 
to the shock arise only from changes in channel area, the Mach number is 
2 function of the area. ‘Taking this over to propagation in the channel 


formed by neighbouring rays, we have the functional dependence 
dA = A(M) (2) 


is the second relation between 4 and .W. ‘This is the only assumption in the 
theory and we can proceed from (1) and (2) without further approximation. 
Qualitatively, the results are independent of the precise choice of (2) provided 
only that 4 is a decreasing function of M. For numerical results we make 
ise of the function A(M) obtained and used recently by R. F. Chisnell 
(1957) for the motion of shocks down converging channels. In an earlier 
paper, Chester (1954) found that for a small change dd in channel area the 
orresponding change in Mach number is given by 


(3) 


vhere A(1) is a slowly varying function, decreasing from 0-5 at VW = 1 
to 03941 (for y 1-4) as W-- x. Chisnell suggested that the integrated 
form of (3) should give a good approximation for a channel of slowly varying 
cross-section; his work on the converging cylindrical shock confirms this 


view. On integration, (3) gives 


2M dM 


) WP KMS “) 


A=kf(M), f(.M) = exp: 
vhere & is an arbitrary constant which may be different for each channel, 
i.e. k = Rk(8). In many cases the shock is initially of constant strength 
with A/ and 4 taking constant values .W, and A4,. ‘Then, for each channel 
k has the same value 4, f(M,) and it can be absorbed into f(./). Even if VW, 
is not constant on x = 0, the k can be suppressed by a careful choice of the 
coordinate 3; it is only necessary to arrange that 4, = f(.M,). For the 
time being we assume this is done. Later, however, we shall see that the 
lependence of k on 8 can arise in a less trivial way. 

The function A(./) is given by 


2 1—p? ! (y —1)M2+2 
F ~ cae We 3 V2 Ry er MER ys. RUM od 
\ (MW) = 2 (1 y4 1 , Je l M } ’ bu 2yM?—(y— 1)’ 


nd its graph is shown in figure 2. Chisnell has shown that the integral 
(+) can be evaluated explicitly; a graph of log,, f(.V/) is given in figure 3, 
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With A = A.M), (1) becomes a second-order hyperbolic equation for VM 
with two independent variables. In the exact formulation of the problen 
there are three independent variables; much of the mathematical simpliti- 
cation is in the reduction of the number of independent variables. ‘The 
solution represents waves moving in each direction on the shock face 





Figure 2. Graph of Chester’s function A(.1/). 


Figure 3. Graph of the function log,, f(.\/) given by equation (4) 


Since equation (1) expresses a kinematic relation, these waves are a furth 


xample of ‘kinematic waves’ in the sense used by Lighthill & Whithan 
(1955). In tact, this example extends the idea of kinematic waves since 


the examples studied previously involved only first-order differential 


quations and, consequently, the wave propagation was in one direction 


only. It turns out that there is a close analogy between the waves travelling 


? 
& 
4 
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m the shock and sound waves of finite amplitude in one-dimensional gas 
dynamics, and all the mathematical methods used in that field are available 
for the present case*. ‘The propagation speed for a wave, i.e. rate of change 
of 8 with respect to x, is an increasing function of M so that waves carrying 
a decrease in the value of M spread out like expansion waves in gas dynamics ; 
similarly, the profile of a wave carrying an increase in M steepens like a 
compression wave. In the latter case, the wave will eventually break and 
then to complete the solution a discontinuity in Mach number and in shock 
slope must be fitted in. This discontinuity is analogous to the shock wave 
of gas dynamics and it must be included in such a way that the appropriate 
‘shock relations’ connecting its speed with the jumps in WV, ete., are 
satisfied. 

Since all the features of gas dynamics arise in the study of the wave 
motion on the shock, it is desirable to use the same terminology because 
it automatically conjures up the right ideas. However, the discontinuous 
wave would then be called a ‘shock’ and unless the word is qualified in 
some way there may be confusion with its direct use for the true shock. 
lo avoid this, we shall always refer to these ‘shocks’ moving on the true 
shock as shock-shocks. 

‘The waves on the shock are interpreted as the trace of cylindrical waves 
which are spreading out in the flow behind the shock. ‘Thus, an ‘ expansion 
wave’ on the shock is the trace of a cylindrical expansion wave in the flow; 
a ‘compression wave’ is interpreted similarly. ‘lherefore, a shock—shock 
must be the trace of a genuine shock produced in the flow behind the main 
shock. ‘Thus, it arises in Mach reflection and represents the three shock 
intersection characteristic of that phenomenon. ‘This feature of the theory 
is particularly valuable since Mach reflection is of great importance in 
diffraction theory and is difficult to treat theoretically in all but the simplest 
cases. All the mathematical details of the wave motion, including a 
discussion of the appropriate conditions relating quantities on the two sides 
of a shock-shock, are given in § 2. 

It is useful to think of the theory as the generalization to shock waves 
of the theory of geometrical acoustics. Geometrical acoustics applies to 
weak sound pulses and is closely related to the special case of weak shocks 
in the present theory. But there are differences, connected with the 
linearization in geometrical acoustics, which can be important. First, in 
geometrical acoustics (1) is replaced by the following determination of 4, 
which is independent of M. Ifthe shock is very weak, its velocity at different 
points varies only slightly and is always close to the sound speed a. In 
accordance with the linear theory of acoustics, these small variations are 
neglected and the propagation speed takes the constant value a). Hence, 
the rays are straight lines, determined once and for all as the normals to a 
ziven initial shock position, and the area A of any ray tube can be calculated. 


* It is assumed in this paper that the reader is familiar with the theory of one- 
limensional gas dynamics as presented for example, in Courant & Freidrichs (1948), 
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he variation of the strength of the shock is then deduced from th 
assumption that the energy flux down any ray tube remains constant along 
the tube. Since the energy flux is proportional to the square of the amplitud: 
of the wave multiplied by 4, it follows from this argument that the strengtl 
is proportional to -'*. For a weak shock, the strength is proportiona' 
to .V/ 1, so that 


We mav note that this agrees with the Chisnell formula as \/ —- 1, sinc 


A = 0-5 in (4). Even for weak shocks, however, geometrical acoustics is 
inadequate in certain cases and the general formulation must be used 
his is due to the linearization which is introduced by assuming that the 
propagation speed can be approximated by ay; although the variations 


the speed are small they cannot always be neglected. Consider, fo 





xample, a shock which is initially concave forward as in figure +. TI 
i 
Ft —— + 
— 
; 
sure 4. Shock positions and ravs according to geometrical acoustics ; AB, C4 


EFG, HIJK represent successive positions of the shock. 


1 
+ 


to the initial surtace form an envelope, called a * caustic’, at whicl 
| 0, and consequently geometrical acoustics predicts infinite strength. 
\loreover, bevond the caustic the position of the shock as calculated by 
seometrical acoustics folds over itself. ‘lo obtain the true behaviour it 1s 
bsolutely essential for the small variations in speed and the corresponding 
listortion of the rays to be included. ‘hen nothing very remarkable 
appens. As the strength increases in the concave part, the shock moves 
taster there and tends to smooth out the shape; the rays are pushed away 
om each other and now avoid any intersection. In fact the region of 
he shock which was originally concave overshoots and becomes convex. 
hen the velocity decreases until the shock ultimately smooths out into ; 
plane. ‘The true picture takes the form sketched in figure 1; the appropriat: 
athematical discussion using equations (1) and (2) is given in $4. The 
ibove argument is also a rough explanation of the observed result that 
plane shocks are very stable; the theory of this paper puts the argument 
in mathematical form. For further discussion of the ditferences betwee: 


eak shocks and sound pulses, reference may be made to the author’s 


paper on weak shock waves (Whitham 1956). 
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lurning now to the question of applications, it is clear from the outset 
that the theory can only be expected to apply to certain types of problem. 
‘he simplicity of the method is achieved by avoiding a detailed discussion 
of the flow behind the shock through the assumption that 4 is a function 
ot MW. Clearly, the flow cannot be forced into such a subsidiary role in 
general. It would not be feasible, for example, in explosion problems in 
which the propagation towards the shock of disturbances originating far 
behind it must be the prime consideration. But, for the diffraction of a 
uniform plane shock by an obstacle, the disturbance originates at the shock 
so that it is more understandable (although still surprising, perhaps) that 
the discussion can be limited to a neighbourhood of the shock. ‘he accuracy 
of the results is a measure of how far this is possible. Again, according to 
the rough argument given in the previous paragraph, stability is largely a 
question of local adjustment of the flow near the shock and can be included 
in the applications. 

Mathematically, the easiest problem to solve is the diffraction of a shock 
moving along a non-uniform wall. If the wall always turns away from the 
How region, an ‘expansion wave’ originating at the wall moves out along 
the shock and the solution is a ‘simple wave’ in the terminology of gas 
dynamics. If the wall turns towards the fluid a ‘compression wave’ 1s 
sent along the shock and eventually ashock—shock is formed. ‘Uhe appropriate 
relations connecting quantities on the two sides of the shock-shock must 
be introduced in this case. A special case is diffraction by a wedge in which 
the solution is given entirely by a shock-shock separating two uniform 
regions; this is the familiar Mach reflection. It must be pointed out 
immediately that the present theory does not throw any further light on 
the difficulties in the conventional solution for conditions at a three shock 
intersection. Its main contribution is to give a method for treating variations 

the three shock configuration which would be caused, for example, by 
further curvature of the wall of the wedge. ‘he special case of a small 
corner (both expansive and compressive) offers an opportunity of checking 
the results with those of the linear theories of Lighthill (1949) and ‘Ting & 
Ludloff (1952). The comparison shows that the approximate method is 

ost suitable for strong shocks with Mach number greater than about 
nd must be used with care for weaker ones. ‘The predicted changes in 
\lach number are in good agreement with Lighthill’s values for all strengths, 
but for very weak shocks the geometry of the disturbed flow is not given 
curately. ‘The details of the application to diffraction problems are set 
in $3. 

In $4, the stability of plane shocks is considered and the results are 
mpared with those found by Blackburn (1953) and Freeman (1955) 
king with Lighthill's linear theory. It is seen that the stability predicted 

by the theory of this paper is achieved by a ditterent proc ss. ‘The essential 
echanism here is the non-linearity of the waves and the dissipation of their 
ie 


y 
s 


v in shock-shocks (in the same way that shock waves dissipate energy 
ordinary sound waves). In the linearized theory of Blackburn and 
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Freeman, the stability depends on damping of the waves through some 
diffusion process (although, of course, this does not appear explicitiy in 
their work). ‘lhe non-linear theory predicts that a sinusoidal disturbanc: 
on a nearly plane shock will decay with time like t-', whereas Blackburn and 
Freeman predict decay as t''? for strong shocks and ¢-*? otherwise. It is 
not clear at this stage which of the two effects is more important or how the 
two sets of results are related. 

Finally, in §5, we continue Butler’s theory (Butler 1955) of the instability 
of a converging cylindrical shock. ‘his work is closely similar to Butler’s 
investigation and the results are the same. ‘lhe advantage of the derivation 
given here is that we linearize the non-linear equation for .V/ by the hodograph 
transformation and no approximation is made. Although Butler does not 
obtain or use (1), his method is essentially equivalent to assuming small 
perturbations about the symmetric solution in equation (1) and retaining 
mly the linear terms in the perturbations. “Vhe method given in $5 is free 
from questions of the validity of the linear approximations and confirms 
Butler’s results. It also offers an opportunity to show how the hodograph 
transformation may be used in this work. 

The relations corresponding to (1) in three-dimensional problems are 
sbtained in Part Il. In particular, they are specialized to axisymmetrical 
problems. ‘Then, as might be expected, the waves on the shock are analogous 
to cylindrical waves in gas dynamics. For example, diffraction of a plane 
shock by a cone is analogous to the problem in gas dynamics of a cylindrical 
piston expanding with uniform velocity; the same method of solution can 
be used. 


(GENERAL THEORY OF THE WAVE MOTION ON THE SHOCK 


‘lo establish the geometrical relationship between J/ and 4, considei 
the curvilinear quadrilateral PORS with vertices (2,8), (%+6z, 9), 
(x+62,8+68), (%,8+68) respectively (see figure 5). Let x, 8) be 
the angle made by the ray with a fixed direction. Since the sides PS 
ind OR are Adp and {A+(eA/ex)dx}58, respectively, and the distance 
between them is ./d5x, the change in ray inclination from P to S is 


OR-PS 1 2A 


—____ = — —— $f, 
PO M ez 
Hence c6 a cd (5) 
cB Max F 


Since the inclination of the S-curves is }7+4, a similar argument shows 


that 
cé P 1 cM 
ae ©" ) 


if # is eliminated from (5) and (6), equation (1) quoted in the Introduction 





is obtained. It is convenient, however, to work with the two first-order 


equations (5) and (6) instead of the second-order equation (1). 








Se aeee 
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We now assume in (5) and (6) that A = A(M) where A'(M) < 0. Then 
20 AM) eM 





Se Rigs — 0) 7 
ae eS (7) 
06 1 oM 

ae = (), 8 
da A(M) op (8) 


As noted in the last section, these equations are like the equations of non-linear 
sound waves and can be treated in the same way. Once the functions 
M(x, 8) and @(, 8) have been found, the coordinates («, 8) may be related 
to the Cartesian coordinates (x,y) through the relations y = f Msin 6 da, 


(A+A,8a)88 





Figure 5. Neighbouring « and f curves in a region of continuous change in M and 0. 


x = [| Mcos@dz obtained by integrating along a ray. ‘The wave property 
is obtained by writing (7) and (8) in characteristic form, in which only 
derivatives in one direction appear. The characteristic form 1s 


3 rl r dM 
ae 7 ae be 
(= sind =)( sae x) ) (9) 


where c is the function of M given by 


_ =a. = a. (10) 


it is easily verified that these equations are completely equivalent to the 
original ones. ‘They show that 


“dM 
6+ | ——- =constant on — =c 11 
J Ac du ‘ ( ) 
i.e. on a wave moving in the direction of increasing B with speed c, and 
*dM dp 
@é— | —— =constant on — = —c, 12 
| Ac dz ; (12) 


i.e. on a wave moving in the direction of decreasing B with speed c. The 
expressions 0+ f{dM/Ac in (11) and (12) correspond to the Riemann 
invariants of gas dynamics. 


F.M. L 
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In the special case of a simple wave moving in the direction of £8 
increasing, 


-dM s 


e= | Ac = constant everywhere ; (13) 
| 2 : 


hence, from (11), @and M must be individually constant on each characteristic 
curve d8/dx = c(M), and this curve is a straight line in the (a, 8) plane. 
‘The solution can then be determined completely in terms of appropriate 
boundary values. For example, let us suppose that @ is a given function 
6,(%) on 8 = 0 and that initially the shock is undisturbed with 6 = 0, 
M=M,. (This example arises in the next section for the motion of a 
shock along a wall.) Using the initial conditions, (13) determines the 
relation between @ and M, and in particular shows that the value /,, of 
M on f = 01s given immediately in terms of 6,, by 


“Mw dM 


° }yy, Ac © 


(14) 
In the («, 8) plane, the slope of a characteristic starting on 8 = 0 at « = «,, 
is c(.M,,(x,,)). Hence, since the characteristics are straight lines on which 
6 and M are constant, we have 


6=6,(%,), M=M,(2,) (15) 
on B = (x—2,,) (M,,). (16) 


Since M,,,is a known function of «,,, equation (16) determines ~,,,as a function 
of x and 8; then (15) gives 6 and M at (z, 8). 

In the general case, however, waves propagate in both directions; then 
(11) and (12) are the basic equations for the well known numerical method 
of characteristics. 

‘To complete the theory we must consider the propagation of shock— 
shocks, i.e. discontinuities of M and @. First, a simple theory is given 
which is suitable when the changes in M and @ are not too large; then, 
we shall see how this should be modified in other cases. Consider the 
neighbourhood of the discontinuity in two successive positions of the 
shock as shown in figure 6 (the rays through the corners are shown as 
broken lines). Let the difference in the « coordinates for the two shock 
positions be Ax and the difference in the 8 coordinates of the rays be Af, 
and let subscripts 0 and 1 be used for values ahead of and behind the 
discontinuity, respectively. ‘Then, in figure 6, PO = M, Aa, OR = A, AB, 
SR = M, Ax, PS = A, Af. Expressing the distance PR in two alternative 
ways, we have 


(A, Ap)? +(M, Aa)? = (M, Ax)? +(4, AB)?. 


But, the ratio Af/A« is the shock-shock velocity C in the (z, 8) coordinates ; 
hence, 
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If it is assumed that the functional relation (4) between A and M applies 
even for the sharp change in channel section at a shock-shock, the velocity C 
is determined by (17) in terms of M, and M,. ‘To deduce the corresponding 
change in 6, we note that angle OPS = 6, — 0, and this angle can be found 
from the geometry in figure 6. For, 


cot(@, — 9) == tan(RPO + RPS) 


_ Ay C+ mS / os. ionh 
M, Ay © M, Ay 


Substituting from (17), we have 





A,M,+4,M 
cot — 8) = Cap waRa(az— Ape ") 
It is easily verified for very weak shock-shocks that the velocity given 
by (17) reduces to the velocity (10) as M,-—> My, A, — Ap, and that, for 
small changes in M and @, (18) gives the same relation as (13). 














Figure 6. Neighbouring « and § curves at a shock—shock, 


Although (17) and (18), together with (4), determine weak or even 
moderate shock-shocks with reasonable accuracy, we must consider further 
the question of stronger ones. We only go into this briefly to give the 
main ideas because considerably more labour would be involved in practical 
applications; this is not worthwhile until it is seen whether the original 
approximations prove satisfactory in practice. Nevertheless, the extensions 
do show up some valuable theoretical points. The limitation on the above 
relations for the shock-shocks is that for sufficiently large jumps in M, the 

L2 
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dependence of A, on M, will not be given accurately by (4), which assumes 
that the channel section varies only slowly. Moreover, it is not just a 
question of establishing the correct formulae relating M and A at an abrupt 
change in channel section. In fact these formulae have been found by 
Laporte (1954). We must remember that there is a third shock in the flow 
behind the main shock, which must be considered in an accurate treatment 
of the conditions across a shock-shock. It does not invalidate the relations 
(17) and (18), but in general we cannot use the simple channel formula for 
the dependence of A, on M,. Of course, it depends on the magnitude of 
the discontinuities; if 1/,—M, is not too large, the third shock will be 
weak and the channel formula can give a reasonable approximation. 

The limitation can be seen in another way. Suppose we accept the 
approximate shock-shock relations so that M,—M, and the velocity C 
are determined in terms of the angle change 6,—6, and the initial Mach 
number M,. Then we can go on to fit in the third shock. But its strength 
must be chosen so that both the pressure and the stream deflection behind 
it are the same as behind the Mach shock. One of these conditions is 
sufficient to determine this third shock and we have a further condition 
still to satisfy. ‘This corresponds to the fact that we have really assumed 
one condition too many in the original shock-shock conditions, namely, 
that the relation between A, and M, is known. If we relax this condition, 
the full three-shock theory will determine the relation in the course of 
fitting in the third shock. ‘Thus, for more accurate conditions we may take 
the relations between M,, A,, 0, and My, Apo, % which are given by the 
conventional theory of the three-shock intersection; in particular, these 
will determine A). 

With this more general determination of the shock-shock, let us go on 
to consider further the question of subsequent disturbances moving along 
the shock. In each ray channel, the variations in A and M are still related 
by (4), but in tracing back to determine the factor k at some known point, 
we can only go back as far as the shock—shock, even if the conditions ahead 
of it are uniform. ‘Thus k = A,/f(M,), where A, and M, are the values 
determined from the three-shock relations. Therefore k is a function of B 
which is to be determined in the course of the solution. Hence, in 
equations (5) and (6) the general form 


A = k(8)f(M) 


must be used. Fortunately, if the new variable [k(8)d8 is taken instead 
of 8, the equations are the same as (7) with A(V) replaced by f(M), and 
they can be solved in a similar way; we omit the details. ‘The whole thing 
is rather like the question of entropy changes in gas dynamics. First of 
all one assumes that the pressure p and density p are functionally related 
(p <p” usually) and this leads to simple waves and so on. But, then, 
since compression waves break, shocks have to be considered and they 
involve entropy changes so that p is no longer a function of p alone; behind 
the shock, the entropy is constant on each particle path. We have the 
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analogous situation with A and M similar to p and p, and k playing a role 
similar to the entropy. Our simple theory of shock-shocks is rather like 
neglecting entropy changes at shocks in ordinary gas dynamics. ‘This is 
known to give quite accurate results if the shock is not too strong, and we 
expect the same to be true here. 

In the applications described in the following sections it will be necessary 
to make use of the properties of the functions c(M) and {(Ac)-'dM. The 
functions are derived from equation (3), i.e. 

dA 2M 
A (M?-1)K(M) 
and graphs of K(M) and A(M) have already been given in figures 2 and 3. 
From (10), we have 
c(M) = (— M/AA’)!? = {3(M?—-1) K(M)}?/4. (19) 
Now c(M) is the propagation speed, i.e. the rate of change of 8 with respect 
to « of a wave in the («, 8) variables; hence the quantity Ac is the rate of 


dM, 


Ac 





2 3 4 5 6 ? 8 9 


Figure 7. Graph of the propagation speed against MM. The full line refers to the 
present theory (equation (19)) ; the broken line refers to the acoustic value 
given by (24). 


change of distance with respect to « since Adf is the line element in the 
direction of 8 increasing. ‘Thus Ac is a more useful as well as more 
convenient quantity to compute; the graph of Ac is shown in figure 7. 
The Riemann variable is given by 

-M dM -M ( in ) 1/2 

| == | \GP=1) Kopf dM, 

J, Ac Jy |(M?-1)K(M)J . 
and its graph is shown in figure 8. 

The approximations of the various functions in the special cases of 

weak and strong shocks are easily obtained using the results that K—> 0-5 
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as M->1, and K > 0-3941 (for y = 1-4) as M+ a. The appropriate 
formulae are 

K(M) ~ 055, Ac ~ (M—1)!2/212, ) 
A (oy M1 | 








Yio, |” 20 ~ 219(M,— 124,’ +} asM—>1, (20) 
“ dM | 
—— ~ 232/(M—1)!2—(M,—1)"23, | 
a ( 2 ( ) | 
and 
K(M) ~ 0-3941, Ac ~ n2M, | 
A ( yg Mri 
—~(—] , c~ no's - -, 
Mu = . MoM as M> 2. (21) 
a. 


MM, Ac _ se M, ; 


where n = 2/K(2) = 5-0743, 
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Figure 8. Graph of the function | dM/Ac. 
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3. DIFFRACTION OF PLANE SHOCKS 


We consider a shock moving along a wall of given shape, specified by 
giving the inclination 6, as a function of the distance s along it. We suppose 
that the wall is straight up to a certain point and that the shock is initially 
moving with constant Mach number V,. The foot of the shock must 
always be normal to the wall; hence, the wall is a ray. If the wall is taken 
as 8 = 0, then @ is given on 8=0. Provided that no shock-shocks occur, 
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the solution to this problem is the simple wave discussed in the previous 
section. Actually @,, is given here as a function of s whereas we require 
its values in terms of « These can be found, however. For, on 8 = 0 





ds 
a= |=, 
J) M,, 
and .W/,,,is determined from 6,, by the simple wave relation (see equation (14)) 
: “Mw dM 
0. = | : 
JM, Ac 


Expansion round a convex corner 

For the special case of a convex corner, 0, jumps from zero to a negative 
value and the solution is a centred simple wave. In the («, 8) plane, the 
characteristics (16) for the disturbed region form a fan, since they all start 
at the same point on 8=(0. The equation of each characteristic is 
8B = ac(M), and since «,. no longer appears this relation determines M 
immediately as a function of 8/x. The shape of the shock will be as shown 
in figure 9; the radial lines shown here correspond to the characteristics 
and on each of them M is constant. The first disturbance spreads out 
on the shock at a rate given by d8/dx = c(M,). If we choose f as the value 
of the distance y from the wall in the initial undisturbed motion so that 
A, = 1, then the actual speed of the first disturbance is a) c(Mo) (since 
a% = at). 








Figure 9. Diffraction of a shock around a convex corner. 


Now for a small bend in the wall, i.e. for @,, small, we can compare 
the results with the linear theory given Lighthill (1949). To be precise 
we compare the values predicted for the Mach number at the wall, /,,, 
and for the speed of propagation of the disturbance. For small @,,, 
(14) reduces to 
M,,,—M, = Ac(M,)6,, = 9,, §3(M2—1)K(M,)}!2, (22) 











160 G. B. Whitham 


(see (19)). We compare this with Lighthill’s value in two extreme cases, 
M,-— 1 and M,-> «. For weak shocks, (22) si 


M,, —M, ~ 6, {3(J M,—- 1)}! 
whereas Lighthill has 8/(37) times this. For strong shocks, (22) gives 
M,,—M, ~ 0:4439M, 84,3 


Lighthill’s value has to be taken from a graph and all we can say is that the 
numerical factor is rather less than 0-5. In view of the relative simplicity 
of the derivation of (22) the results are remarkably good. 

‘Turning now to a comparison of the speeds of propagation, it follows 
from the theory of sound that the first possible disturbance travels out in 
the flow behind the shock with the local sound speed a relative to the flow 
velocity u. ‘Therefore the disturbance travels along the shock with speed 

(a?—(U—u)?}}2, (23) 
where U is the undisturbed shock velocity. The quantities U, a, u can 
all be expressed in terms of M, and it is found that (23) is a,c*, where 

2 Apacs 2 \ 2 

a — | R= = 1) MR + 2p (24) 
(+ 1)Mi 
This is to be compared with the speed cy = c(M)) given by (19). Since 
Ay = 1, ¢ is the same as the graph of Ac in figure 7; the graph of c* is also 
shown in figure 7. For weak shocks, 
Cy (4(M, —1)})", CF me {2(M, — 1)}"; (25) 


for strong shocks 





Cy ™ 0-4439M,, c* = 0-4083M). (26) 


Thus, the values are in reasonable agreement for shocks with M, > 2, say; 
but, for very weak shocks cy = $c*. ‘This discrepancy arises because the 
present theory cannot avoid concentrating the change in M over a relatively 
small part of the shock. For the stronger shocks, Lighthill’s work shows 
that this concentration is correct; in fact he finds that the curvature becomes 
infinite as /,-> oo. But for weak shocks the disturbance should be spread 
out over the whole of that part of the shock which is inside the sonic circle. 
The present theory compromises by concentrating the disturbance halfway 
to the sonic circle. Generally speaking then, the theory is more suitable 
for moderately strong shocks. However, the correct prediction of the 
value of M,—,, even in the weak case should not be ignored. 

When the magnitude of 6, is not small, MW, must be found from the 
exact form of (14) using the graph shown in figure 8. The solution is the 
centred simple wave determined from c(M) = 8/x, c(M,,) < B/x <¢c(M)). 
As an example, the details of the solution are given for an infinitely strong 
shock. First of all, from (14) and (21), the Mach number at the wall is 
given by 


M,, = Myexp(@,,/vn). (27) 
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Then the region c(M,,) < B/« <c(Mp)) covered by the simple wave can be 
written 
exp{(m+1)),,/vn} < Bvn/aM, < 1. 
In the simple wave (using the results in (21)) 
Ages 


_f. es M, 
c(M) = ame a and 6 = —vnlog ve 


Therefore, 


M - e : | 
M, \aMy | | 
. ! » exp{(nt+1)8,/vn}< ah <4, (28) 





| in | aM 
a coe log hie . 
n+1 aM, J 
Along the shock, ex/AdB = —sin@, dy/AdB = cos@; therefore, at time 
t = a/a», the shock is given in terms of the parameter 8 by 
B 
x = «aM,,cos@,,— [ (M,/M)" siné dp, 
ie (29) 
y=aM,,sin8@,,+ | (M,/M)" cos @ dp. 
~ 0 


It should be noted that x/M,« and y/M),« are functions of the single 
quantity 8/M,«, so that the shock pattern expands uniformly with time and 
a change in M, involves only a change of scale. The values of x and y in 
the simple wave are most easily calculated from (29) with @ as the parameter 
instead of 8; they are 





x (n+1)%? 4, -: } 
= —_,— e9V" sin(A— 9), | 
Mya ee Ne | 
eee 0,<0<0 (30) 
y a+ly* ... 
Mya = —— ey cos(A — 8), | 
where tanA = vm. The shape of the shock is plotted in figure 10 for the 
special case 0,, = — 47. 


For weaker shocks there is a limit on the magnitude of 6,,.. For, M,, 
cannot decrease below unity; hence, if @,, decreases below the value 
Gm given by 

9 (" dM , 


lim = , 
m, Ac 


the solution breaks down. Presumably this corresponds to separation of 
the flow at the corner which is known to occur in certain cases. 


Compression at a concave corner: diffraction by a wedge 

For a concave corner, the solution given by this theory is a shock-shock 
separating two regions in which M and @ are constant (see figure 11); this 
is Mach reflection. Following the remarks in §2, the most accurate 
determination of the solution in this theory is the conventional three-shock 
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intersection treatment. ‘Therefore, the only question that need be discussed 
here is the extent to which the conventional treatment can be replaced by 
our weak shock-shock relations (17) and (18). We consider the values 
predicted by (17) and (18) for the angle, y, between the wall and the line 
joining the corner to the triple point. When the triple point is at (a, £), 








0.27 








ITTF, 


a a 


a 





’ 
uo 





‘ 
o 
Sie (eG are 





Figure 10. Diffraction of a strong shock around a right-angled corner ; shape of 
shock calculated from equation (30). 








Figure 11. Diffraction of a shock by a wedge. 
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uf 


tan y = M,,8/A,,%, and B/x is equal to the shock-shock velocity C. 


Therefore, from (17), 





A,, {1-—(M,/M,,)?) ¥? 
tany = — <———— > > 31 
" Ay | 1-(A,,/Ao)? J GN) 
Also, from (18), 
M, 1+A,,M,,/A,M, (32) 


Go = ass , 
corte = Mf, (1 —M3M2\(1 — Az Aye 


For strong shocks, the curve of y against 0,, becomes independent of Mp, 
since A,,/Ay is a function of M,,/M, (equation (26)). This curve is drawn 
in figure 12 and compared with the corresponding curve obtained from 
the three-shock theory (assuming that the Mach shock is approximately 
straight). We expect our shock-shock relations to apply when @,, is small 


Figure 12. Values of y and @,, for the diffraction of a strong shock by a wedge. 
The full line refers to the present theory ; the broken line refers to the three- 


shock theory. 


but to diverge ultimately from the more accurate value. It is seen that the 
error for small 6@,, is about as expected, being of the same order as the 
discrepancy in (24). Then, fortuitously, the curves come closer together 
and actually cross. However, in the three-shock theory there is an upper 
limit on @,, at which Mach reflection goes over into regular reflection, while 
the simple shock-shock relations become useless as they continue to predict 
Mach reflection. It is perhaps worth noting that if A,, is related to M,, 
by Laporte’s formulae for a finite change in channel section, a cut-off value 
for @,, is obtained (the value being @,. = 33-6°). But as explained in § 2, 
this is not the correct direction in which to look for more accurate shock—shock 


conditions. 
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As an example of a shock of moderate strength the graph of y against 6,, 
was found for M, = 2-42, corresponding to a pressure ratio 20/3. This 
is shown in figure 13 and compared there with the experimental values 
quoted in Bleakney & Taub (1949). 


° ae ne .° ane 
0° C 2c 3¢ 40 5 


Figure 13. Comparison of the graph of x against 6,, with the experimental results of 
Bleakney and Taub (indicated by broken line) in the case My 2°42. 


Wall of arbitrary shape 

For the general problem, the solution must be obtained as in the 
analogous problems of gas dynamics; the slope of the wall corresponds 
to the velocity of the ‘piston’ in gas dynamics. If shock-shocks are not 
formed the solution is the simple wave described by (15) and (16) of § 2. 
Weak shock-shocks can be fitted in by a technique developed for gas 
dynamics and supersonic flow. ‘They are formed when the characteristics 
overlap leading to a multivalued solution in (15), (16), because for points 
(x, 8) in the overlapping region there will be more than one value of ~,,. 
This is avoided by fitting in a shock-shock, according to the relations (17), 
(19). "The rule (Whitham 1952) for fitting in a weak shock-shock is obtained 
in terms of the values of «,, on pairs of characteristics which meet at the 
shock~shock. Let the two values of «,, for such a pair of characteristics 
be «, and x. ‘Then, the values of x, and x, are calculated for each 8 from 


the two relations 


a: a 33 
F(a) — F(a) sa 
|” F(x) dx = 3(x.—%){F (2) + Flaa)} (34) 


where F(x) is the function 


1/¢g — 1/c(M,,()). 
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For the derivation of this type of result and for detailed explanations of the 
method of using it see Whitham (1952) (especially the Appendix). It is 
consistent with the approximations made already to approximate F(«) by 
the expressions 
c(M,,(«)) —¢ dc 

F(a) = as == {M,,(«) — My}co? law |. al (35) 
At each line 8 = constant in the («, 8) plane, all the characteristics with 
their starting points «,, in the range «, <«,, <a, are omitted since they 
have already been cut off by the shock-shock. In this way all except one 
of the values «,, given by (16) will be excluded, and the solution becomes 
single-valued. At the shock-shock, F jumps from the value F(«,) to the 
value F(«,); the corresponding jump in M is then given by (35). Certain 
applications of these results will be made in the next section. 


Motion of a shock between two walls 

The motion of a shock between two walls may also be treated by this 
theory. It would be analogous to the problem in gas dynamics of the waves 
produced in a tube of finite length by pistons in each end of the tube. The 
given shape of the walls corresponds to given motions of the pistons. In 
this case waves would move in each direction across the shock face, being 
reflected from each wall in turn. ‘The solution would require numerical 
computations (by the method of characteristics, say). 


4. STABILITY OF PLANE SHOCKS 

It is well known that plane shocks are stable; that is, if the shock is 
disturbed slightly from the plane shape the disturbance decays as the shock 
propagates. This property is now investigated using the theory of wave 
propagation on the shock, and the rate at which a disturbance will die out 
is determined. 

As a special case we may consider the disturbances generated as the 
shock moves along a wall, and the results of the last section may be used 
to study the decay of the disturbance. Let us take first the case of a small 
bump on the wall. Then @,,=0 except on a finite length of the wall. 
In this case it can be shown that after a sufficient time the disturbance will 
become an ‘ N-wave’, i.e. there is a shock—shock at both the head and tail 
of the disturbance. The maximum values of 6 and M— M, are taken at 
the shock-shock, so it is sufficient to calculate the rate at which these values 
decay. Since the wave is moving into an undisturbed part of the shock, 


F(«,) = 0 in (33) and (34), and therefore 


oe ab. ee a 
p= F(a) ’ | ; F(a) dx = 3(a, a1) F(x.). 
On eliminating (%.—.,), we have the relation 
2 


= gates | Blaha 


B 
(a2) 0 
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to determine the characteristic variable x, behind the shock in terms of 8. 
For large f it is clear that the corresponding x, tends to the zero of F. If 
this zero is denoted by % we have 

“Xe 1/2 

F(a) ~ 45 | F(x)dx> . (36) 


(PR 40 J 


The changes in M and @ at the shock-shock are proportional to F(x,): 


from (35), 
M-M, = E = | F(a), 


and, trom (13) with A, = 1, we have 


= (M-— M,)/¢e. 


Thus, the disturbance decays like 8~' as it moves away from the wall, 
and since the disturbance is weak, 8 is approximately equal to the distance 
from the wall. 

Blackburn (1953) investigated the case of a shock moving along a 
sinusoidal wall using Lighthill’s linear theory. According to the present 
theory the sinusoidal variation sends out a series of successive compression 
and expansion waves, and the compression waves eventually break to form 
a series of shock-shocks. As 8 > ©, the values of x, and a, given in (33) 
approach successive zeros of F(x), say % and %)+/, and we see from (33) 
that the jump in F(z) across each shock is given by 

AF = —//B. 
In place of /, we may introduce the wavelength A of the sinusoidal wall 
which is approximately IM,; also 8 = y. ‘Then the changes in M and 0 
are given in terms of AF as before, and we have 


Cc” dl VU A C dM A 
M = £0 ee fee —, 37 
onelim My z dc Lo = M, y’ ” M | dc iq - MY oa 


It is interesting to note that these values are independent of the amplitude 
of the sinusoidal variations in the wall slope. ‘The factors multiplying A/y 
in (37) are increasing functions of M, so that the stability decreases with 
increasing Mach number. 

Now the decay like A/y predicted here does not agree with Blackburn’s 
results. He finds that the decay is proportional to (A/y)*? in general, but 
as M -» o the law changes over to decay proportional to (A/y)'*. But the 
mechanisms of decay in the two theories are completely different. In the 
present theory, non-linearity is essential; the decay is brought about by 
the shock-shocks. Indeed if the theory were linearized in the usual way, 
the waves would not decay at all. On the other hand, Blackburn’s theory 
is a linear one and he finds a decay due to a damping of the waves 
which does not appear here. If we pursue the analogy to gas dynamics, the 
present theory corresponds to the non-linear theory of sound waves, 
neglecting viscosity and heat conduction except where their effects are 
concentrated in shocks; the Blackburn theory corresponds to the linear 
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theory of acoustics with viscosity and heat conduction included. However, 
it would be dangerous to push this analogy too far, and it should not be used 
in assessing the relative importance of the two mechanisms of stability. 

Another formulation of the stability problem which leads to similar 
results is the initial value problem in which the shape and Mach number 
of the shock are given at some time, i.e. @ and MV are prescribed functions 
of 8 on x=0. ‘The problem of a non-plane piston moving forward with 
uniform velocity into a gas at rest is a special case of this. In this case 
waves move in each direction on the shock. If the shock is initially plane 
and uniform except in a section of finite length, there will be an interaction 
region at first, but eventually the disturbance will separate into two simple 
waves one moving in each direction. Each of these simple waves is similar 
to the solution obtained for the problem of a bump on the wall, and decays 
like ¢??. An accurate solution in the interaction region could be obtained 
numerically using the method of characteristics. But if the disturbance 
is not too large and the interaction region not too wide, the linearized form 
of the solution will give quite good results and it can be given explicitly. 
Let us, for example, consider the case in which M takes the constant value 
M, on « = 0 and 6 is a given function of 8. In the linearized theory, (7) 
and (8) are approximated by 








fale) 1 oM 

a* aia om 
06 ol 
oS ee el (39) 


dn yi Op 
Then, 6 and M both satisfy the wave equation and we have the solution 
0 = g,(B—cy%) +¥2(B + cy %), 
where g, and g, are arbitrary functions. (This solution corresponds to 
approximating the characteristics in the exact theory by lines B + cy«.) 
Since M is a constant on « = 0, (39) shows that 06/¢x = 0. Therefore 
the functions g, and g, are to be determined from initial conditions 6 = 4,(8), 
say, and c#/éx = 0 ona=0. These determine the solution as 
9 = 319 9(B — cy) + Oo(B + co x)}. (40) 
The curvature of the shock is proportional to 26/08 and so satisfies a similar 
rule. ‘This shows very directly the tendency of the curvature of the shock 
to be averaged out. For example, the curvature never exceeds the maximum 
value of the initial curvature. Thus even the linear theory rectifies the failure 
of geometrical acoustics near a caustic, and the curvature remains bounded. 
The linear theory is not uniformly valid as « becomes large due to the 
divergence of the characteristics, and we must go over to the accurate 
determination. 

Freeman (1955) considers the special case of a shock which is initially 
sinusoidal in shape, and as would be expected his results for the decay of 
the disturbance are similar to Blackburn’s. In our theory we have a rather 
complicated situation with shock-shocks moving in each direction and 
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interacting with each other. But the results will be essentially the same as 
(37) with the disturbance decaying like 1/t. 

As a preliminary investigation Freeman considers the shock produced 
by a nearly plane wedge-shaped piston moving forward with uniform 
velocity. It is perhaps worth noting our solution for this case. The initial 
Mach number M, of the shock is constant, and is easily found in terms of 
the piston velocity from the shock conditions. ‘The initial values of 6 are 
6, = —din 8 < 0, 6, = +4din 8 > O, where z— 26 is the angle of the wedge. 
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Figure 14. Shock patterns produced by wedge pistons. 


If the wedge is convex forward, the solution is given by two centred simple 
Through the simple wave moving in the direction of £ increasing 


waves. 
‘i ‘M dM 
weet he 
Ju, 44¢ 
and through the other wave 
: ‘M dM 
G+o0= —- | es, 
m, de 


In between the simple waves, both these relations hold; hence, 6 = 0 
and the Mach number .W, is determined by 


Mi dM | M-M, (41) 


y 
20o0o=> = —_ = 
/ My Ac Agcy 


If the wedge is concave forward, 6 is negative and the solution is given by 
two shock-shocks. ‘The physical plane and the (x, 8) plane are shown for 


each case in figure 14. 
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Although these solutions have been introduced for a rather artificial 
problem, they also furnish results for the more interesting problem of 
reflection of a plane shock from a nearly plane wall. It is only necessary 
to choose the frame of reference so that the air behind the incident shock 
is at rest. Then the reflected shock is given by the above solution; the 
flow pattern for both a convex and a concave corner are shown in figure 15. 








Figure 15. Reflection of shocks at wedge-shaped walls. 


The appropriate values of 5, Mp), a) are obtained from the Mach number M’ 
of the incident shock, the wedge angle 7—28’, and the sound speed a’ 
ahead of the incident shock, using the well-known results for regular 
reflection. When the incident shock moves with Mach number M’ into 
air at rest, the particle velocity wu’ behind it is given by 


i : 1 
@ y+ (m se ir) 


The velocity —w’ is superimposed on the flow in order to obtain the required 
frame of reference. In this frame, the reflected shock moves into air at 
rest and it can be described by the above theory. The sound speed ay is 
given by ‘ 
a {2yM—(y—1)}{(y—1)M? +2} 


a (y+ 1)?M’3 A ’ 
and the Mach number M, and the inclination 6 of the ‘undisturbed part’ 
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of the reflected shock are given by 
{37-1 os 8, a in 
ly+1 — CEST Re a \lytl (yt+1)M" 


5. STABILITY OF A CONVERGING CYLINDRICAL SHOCK 











In contrast to plane shocks, it is found that converging cylindrical and 
spherical shocks are unstable. As noted in § 1, a theoretical demonstration 
of this instability has been given by Butler (1955). We now give a somewhat 
improved presentation of Butler’s theory which avoids making the small 
perturbation approximation. 

A converging cylindrical shock will ultimately become strong, so that 
for the stability investigation it is sufficient to consider the case of strong 
shocks. We suppose that the Mach number is constant and equal to M, 
at some initial time, and we choose f to be the distance of the ray along the 
initial position of the shock so that Ay = 1. ‘Then we take 

A =(M,/M)", n = 5:0743, 
as given in (26), and the equations (7) and (8) for 6 and M become 
06 a nM" oM 


in 2 

ap t Met oa = (#2) 

oO a0 M" 0M 2 

iM 2B = 0. (43) 

For a shock with iitibelas symmetry and initial radius Ro, 6 = — B/Ry 
and M is a function of «. From (42), we see that 

- (2+ 1)Mga\-m+) 
m-m,(- Steep", (44) 


« being chosen so that it reaches the value zero when the shock gets to the 
centre. This is Chisnell’s approximate form of the exact Guderley solution 
(Guderley 1942). For y = 1-4, (n+1)7 = 0-16463; the corresponding 
exponent in Guderley’s solution (calculated with great accuracy by Butler 
(1954)) is 0-16478. In his discussion of stability, Butler chooses the value 
of n from Guderley’s solution. 

Now the question is whether small deviations in the initial shape of the 
shock increase or decrease as the shock contracts. ‘The hodograph trans- 
formation will be used to investigate this question. First it is convenient 
to introduce new variables : 

= (M/M,)**}, x = (n+1)0/vn, s= Mya/vn. 
Then (42) and (43) take the neater form 


oq 2 . . 
+ - 0 
F +a 5g =O (45) 
cx 24 
a, + ap = 0. (46) 


B. 


The symmetrical solution is the one in which q < 1/s, x x 





(It may be 
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noted in passing that another exact solution is 
q=B/s, x = log(B/s); 
this describes the motion of a shock down a certain curved channel.) In 
the hodograph transformation, the roles of the dependent and independent 
variables are interchanged; the transformation formulae for the derivatives 
Fr a= Js, ) —JB,, te >= —Js,, qs = JB, 
where the Jacobian J = oq, x)/e(s, 8). On substituting these expressions 
in (45) and (46), we have the /imear equations 
B,+4?s, = 9, B, +s, = 0. 
It is convenient to eliminate 8 to get the single equation 
Be ee es 
T 50g T Bq = Sy 
Solutions to this equation are 


= gte'", 
where p= —$F(t—m?)?. 
If m = 0, » = —1 gives the symmetrical solution. If m > 1, u is complex 
with A(uw) = —}. ‘Therefore, when g-—> © (as the shock contracts to the 


centre), the harmonics eventually dominate the symmetrical mode. Hence 
the shock is unstable. 

‘The appearance of an imaginary part in yp is also of some interest since 
it shows that the disturbance again takes the form of waves travelling round 
the shock. When the disturbance grows large, it is possible for the Jacobian 
J to vanish. This means that the mapping from the hodograph plane (q, x) 
to the (s, 8) plane ceases to be single-valued, and so-called limit lines appear. 
These correspond to the formation of shock-shocks. When this stage is 
reached, further calculations of the motion would be more easily carried 
out directly in the (s, 8) plane by numerical methods. 


REFERENCES 
BLACKBURN, D. 1953 Ph.D. Thesis, University of Manchester. 
BLEAKNEY, W. & Taus, A. H. 1949 Rev. Mod. Phys. 21, 584-605. 
Butter, D. S. 1954 Armament Research and Development Establishment, Rep. 54/54. 


Butter, D. S. 1955 Symposium on Blast and Shock Waves, Armament Research 
and Development Establishment. 

Cuester, W. 1954 Phil. Mag. (7) 45, 1293-1301. 

CHISNELL, R. F. 1957 3. Fluid Mech. 2 (in the press). 

Courant, R. & Friepricus, K. O. 1948 Supersonic Flow and Shock Waves. New 
York: Interscience. 

FREEMAN, N. C. 1955 Proc. Roy. Soc. A, 228, 341-362. 

GuperLey, G. 1942 Luftfarhtforschung 19, 302-312. 

Laporte, O. 1954 Los Alamos Scientific Laboratory, Rep. no. LA-1740. 

LIGHTHILL, M. J. 1949 Proc. Roy. Soc. A, 198, 454-470. 

LIGHTHILL, M. J. & WuitHaM, G. B. 1955 Proc. Roy. Soc. A, 229, 281-317. 

Tic, L. & LupLorr, H. F. 1952 ¥. Aero. Sci. 19, 317-328. 

Wuitnuam, G. B. 1952 Comm. Pure Appl. Math. 6, 301-348. 

Wuitnam, G. B. 1956 ¥. Fluid Mech. 1, 290-318. 








Damping of surface waves in an incompressible liquid 
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SUMMARY 

The damping of surface waves of small amplitude in liquid 
contained in cylinders has been calculated. Viscous dissipation 
in an assumed laminar boundary layer was taken to be the primary 
cause of damping. Experimental results were obtained for the 
logarithmic decrement as a function of the ratio of liquid height 
to cylinder radius for several water-filled cylinders. Theory 
and experiment were found to be in good agreement. 


I. INTRODUCTION 

The calculation of the natural frequencies of the surface waves in a 
liquid contained within solid boundaries except for one free surface is a 
well-known and largely solved classical problem (see Lamb (1945), § 257); 
but the corresponding problem of the damping of these waves does not 
appear to have been so thoroughly treated. Some early work on the 
viscous damping of surface waves is reported in Lamb (1945, §§ 348 & 349), 
but this was mainly confined to progressive waves far removed from solid 
boundaries. ‘The first attempt to account for the effects of solid boundaries 
on wave damping appears to have been made by Boussinesq as early as 
1878. His theory extends to both progressive and standing waves, and 
was used in more recent times by Keulegan (1948) as a starting point for a 
calculation of the attenuation of solitary waves. Other modern work on 
the damping of progressive waves includes that of Biesel (1949), which 
deals with waves in a channel of finite depth but infinite width, and that of 
Ursell (1952), which concerns the dissipation in the vicinity of vertical 
walls when the depth is infinite. Furthermore, Hunt (1952) has calculated 
the combined effects of finite width and finite depth. The latter authors 
all employed boundary-layer approximations, which are applied in this 
paper to the case of standing waves. In particular, we shall give the results 
of calculations of the damping of standing waves in right circular cylinders. 

Since the calculations involve a number of idealizations, it seemed 
advisable to check the theory with experiment. Accordingly, the damping 
was measured for water in cylinders of different radii as a function of the 
ratio of water depth to cylinder radius. ‘The apparatus used and the 
results obtained are described in $III. In $IV it is shown that, subject 
to certain limitations, the agreement between the theory and experiments 
is satisfactory. 
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II. ANALYsIS 

In the interests of simplicity we have restricted ourselves to the considera- 
tion of small amplitude oscillations. In addition to linearizing the equations, 
this restriction frees us from worries about turbulence. Rough calculations 
indicate that with cylinders of the radii used in the measurements the flow 
is always far from the turbulent state. Some idea of the limits of validity 
of the small ampiitude approximation can be obtained from the results of 
the experiments (§III) by noting what amplitudes of excitation give rise 
to a simple exponential decay law. 

For completeness and to provide a basis for the later work, we first 
sketch the calculation of the natural frequencies. 

Consider a rigid right circular cylinder of radius R with base at s = — $h. 
The equilibrium free surface of the liquid in the cylinder is at z = +3h. 
For liquids with low viscosity, such as water, we can expect that a good 
zeroth approximation is obtained by neglecting viscosity entirely and 
describing the liquid velocity distribution by means of the velocity potential 
function ¢,. Thus, we take 


u= —V¢., (1) 
V7, = 9, (2) 

and Ch, 
— =e (3) 


on rigid boundaries, where n is the normal vector. 
In cylindrical coordinates (r, 6,2) solutions of Laplace’s equation 
subject to the boundary conditions (3) are 


(r,6)\ 
— ne 1 Am 
$,= feoshhy,(=-+ ml}, (4) 
where co aa Scos s6) " 
1 ria = Nns sin s0 foes?) (5) 
Here, r _ Jy pez __ 12 
N ne tie Jo ms R\(1 i m)} (6) 


is chosen so that the integrals over the cylinder cross-section of y?, and 
9 . rr - 
y, are unity. The k,,, are defined by 


J'(Ry,¢R) = 0. (7) 


ms 
Expanding ¢, in terms of these proper functions, we have 
dp, = > {cosh Ring(2 5 bh (Aa, Xms 7 Mas Pins} (8) 
If we expand 7(r, 6) (the height of the free surface at (7,4) above the 
equilibrium plane) in terms of the y,,, and y¥,,, we have 
) 1 
n(r, 0) = > {Ime Xms +Pins Pst (9) 


ms 
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The free surface condition 
=~ =u, (10) 


permits us to express the expansion coefficients 4,,,, B,,, in terms of the 


ms 


‘surface coordinates’ g, ., in the form 
ms 


ms 


r 


(A, Ps = Ins | r 
a" “yy | Hey sinh Ry he (11) 


In these coordinates the kinetic energy becomes 


ms) ms) 


Pes ' Coes -« Conk hk 
T=| tpu?dV = 3p > [¢,.+bi,s5] ——— > (12a) 
MF i me Ras 
where V and p are the liquid volume and density respectively. 
The potential energy 2 of the liquid is the sum of the gravitational 
potential energy ({2,) and the surface tension energy ({2,), where 





OQ, =pg| 2dV = spe ¥ [dine t+Pruah (13 a) 
. v ms 
j on \2 On\2 
and 0, [sof (Py (22) 
JSt- OX oy 7 


bo YR Ge + Pisl: (14) 


(Here S, is the equilibrium free surface z = $A and a is the surface tension 
coefficient.) 

The equations for the normal coordinates q,,, OF P,,, obtained by 
substituting the energy expressions (12a), (13 a) and (14) in Lagrange’s 
equations (see Lamb (1945), $135), are 


pq,(cothk,, .A)/k,,. + (pg 4 ok? .)ding —4 0 F (15) 
From this there result the proper frequencies 
prop | 
w, = gk, .tanhk,, ,Afl + ok®,,/(pg)}. (16) 


The contribution of surface tension is unimportant for the lower modes of 
vibration (k,,, R small) except in cylinders of quite small radius. Since it 
adds only 1°, to the frequency of the mode considered in the smallest 
cylinder used in the measurements (R= 1-5in.), it will henceforth be 
omitted, and Q will be approximated by ©, (given by (13 a)). 

In the above development there is, of course, no damping. To find 
this we turn to the linearized Navier-Stokes equation for an incompressible 


fluid 
in = — ¥(¢e + ’) +vV?u, (17) 
ct i p 


with the equation of continuity 
V.u=0, (18) 
and with the boundary condition 
u=0 (19) 
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on fixed boundaries. Here p is the pressure and v = p/p the kinematic 
viscosity. 
Let us represent the velocity vector in the form 
= —V¢+VxA, (20) 


where ¢ and A are respectively scalar and vector functions of (r, 6, 3, Z). 
It is readily demonstrated that this, together with 


: +2 = = +const., (21) 
constitutes a solution of (17) and (18), provided that 
V7¢ = 0, (22) 
oA 
vwW7A = i (23) 


If the boundary conditions of the problem can be satisfied by means of 
the expression (29) subject to (22) and (23), we may infer that this is the 
complete solution. ‘The boundary conditions at the rigid surfaces S 


(-—V¢+Vx A), = 9, (24) 


together with (22), (23) and the surface condition (10) make possible, 1 
principle, an explicit, rigorous solution for ¢ and A in terms of the surface 
coordinates. However, the general experience gained with boundary 
layers in the last fifty years suggests that this is hardly necessary. It is to 
be expected that ¢ is essentially the velocity potential ¢, characteristic of 
inviscid flow (i.e. the potential for which (¢¢,/én), = 0). 

Let us write 6 = ¢,+¢;. We shall see that the contribution of the 
additional part ¢/ to the velocity is quite small. Assuming this (subject 
to later verification), we have as approximate boundary conditions on the 
rigid surfaces 

[n x (V x A)], = (nx be = (25 a) 
and (n.Vx A), = (0¢;/¢ (25 b) 


Since ¢, is now assumed known, (23) and pe a) determine A. Having 
obtained A, (22) and (25 b) determine ¢/. The method of calculation is, 
in principle at least, one of successive approximation. Starting with the 
zeroth approximation (A = 0, ¢ = ¢,), we calculate a first approximation 
to A. With this a correction to ¢ can be found. In principle this process 
could be repeated to an arbitrary accuracy. In practice, however, the 
first approximation is sufficiently accurate by itself. 

The solution of the equations for A is enormously simplified by noting 
that the region of non-vanishing A is confined to-the immediate vicinity 
of the boundaries. Indeed, assuming the special case of a simple harmonic 
time dependence (A ox e~“”’), we obtain 


(V2+7/P)A = (26) 
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where the boundary-layer thickness / = 1/(v/w) is of the order of 0-1mm 
for the values of v and w arising in the experiments described later. Outside 
this layer, A vanishes exponentially. Consequently, to find A near a 
given boundary, the other boundaries can be ignored. For a standing 
wave this case represents a forced oscillation, requiring a supply of energy 
to counteract viscous dissipation. However, for free oscillations in which 
the damping is small, the frequency w is large in comparison with « in the 
exponential decay factor e~“; thus we may again assume a boundary layer 
of thickness 4/(v/w), and may make use of the suggested approximations. 
In the body of the liquid, including the vicinity of the free surface, A can 
be taken to be zero in this approximation; accordingly, the boundary 
conditions at the free surface are satisfied in terms of ¢, alone. 

To compute the viscous dissipation we use the well-known result 
(see Lamb (1945), § 329) 

dE/dt = —2F 27) 

where E is the total energy (7+), and 2F is the dissipation function 
given by 


2F=p| (VxudVep| (m.Vut)dS—2n| n.ux(Vxu)dS, (28) 
a J §’ J Ss’ 


s 
in which the two surface integrals extend over the whole boundary 
S’ = S+5S, of the liquid volume V. ‘The bars denote averages over a 
cycle for an assumed harmonic time dependence. Equation (28) simplifies 
on noting that u=0 on rigid boundaries, and V x u = 0 approximately 
at the free surface. ‘Therefore 


2F = 2F,+2F, (29) 
approximately, where 
2F = p | : (VxVxA) dV, (30) 
me 2F,--p( (v4)? 4S, (31) 
J Son 


where A is expected to be insignificant away from the solid boundaries. 
The latter integral, which is a simple transformation of a volume integral, 
represents the dissipation in the body of the liquid. It may at first seem 
inconsistent to retain this term, which is clearly of smaller order of magnitude 
than the first term under the assumptions of the present approximation, 
while neglecting the dissipation, due to rotational motion at the free surface, 
represented by the third integral in (28). (This integral measures the 
additional dissipation in the ‘boundary layer’ which must exist near the 
free surface in order that viscous stresses are balanced at this surface— 
where the only possible non-uniform stress is that of surface tension.) 
However, the dissipation in the boundary layer at the free surface is in 
fact of lesser order of magnitude than that in the body of the liquid (see 
Ursell (1952), p.94), and the term (31) gives a reliable indication of the 
magnitude of the dissipation occurring away from the solid boundaries. 
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To find the damping of a single mode, we keep only one term in (8). 
The kinetic and potential energies are then 


‘ } 
, cothk,, A 


T= bi, (12 b) 
Q = bpegins: (13 b) 

For damped oscillations we have 
Gms = Tne emt. (32) 


Assuming « < w (which is justifiable a posteriori) we obtain 


dE = 





Ima (— xe. (33) 


Inserting ¢, from (8) into (31) we obtain after some elementary integrations 
—— 2 #2 —2at 
2F, st 2uRins PE4ms —s (34) 


The differential equation for A in the vicinity of a given wall can be inte- 
grated readily in terms of known functions (assuming the other boundaries 
to be at infinity). However, even this is not necessary. ‘The boundary 
layer thickness / is small compared to the radius of the cylinders considered. 
Hence we can neglect the effects of curvature. (The error made is then 
of the order of //R, which is less than 1%.) With this approximation the 
vector potential in the vicinity of the side walls is 


A, = D,e-v"**" cos sO sinh k,,,.(2 + 3A), 








A, = Dye~V¥“*®” sin s6 sinh R,,,,(2 + 4h), (35) 
pas | 
where ae. —— ] 
** kn R TAs R)( — sinh k,,,. i) 
l ae fa | 2 (36) 
Dy = sphae The) gas) 1a ea)> | 
T= y/(v/ay,5)- | 


With this expression for A it can be seen from (25 b) that the resulting 
¢. is O(lk,,,) = O(l/R). ‘This then verifies the earlier statement that the 
contribution of ¢, to the fluid velocities is small. Indeed, we see that the 
approximation method adopted is essentially an expansion in powers of //R. 

Correct to first order in //R, we find from (35) in conjunction with (30) 
the contribution (2F,) of the side walls to the dissipation due to the rigid 
boundaries to be 
po, ete 1 


~ kk [1—(s/k,,.R)2}] 2V21 


ms 





9 | h 9 
x {(coth ys LL + (she RY = ( FE I~ lh RI. (37) 


ms 





178 K. M. Case and W. C. Parkinson 


The contributions to the dissipation from the tank bottom are obtained 





similarly. ‘Thus, for the vector potential in the vicinity of z = — 3h we 
find 
, Ss ) 
A. = Ee-Vitet+ ton J (Ring 7) COS SO, 
r k,,.R A ms ) | 
1\d ; ‘ (38) 
A, = —Ee-viet" | —_ |_ J (k,,,r) sins, | : 
2 - dr a( nts )s 
A, «®, } 
rn —1 f —Gm,N 
where o. = t- ms *Nins\ (39) 
vi \sinhk,, A 


Calling the contribution to the dissipation 2F, we have (again to first 
order in //R) 
e k #2 gp—2at 
2F ae be ms Ims ve 
a es me "7 ea ca.) oe a 
v2/sinh 2k,,.A 


ms 


(40) 


Inserting the results expressed in (33), (34), (37), and (40) into (27) 
yields an equation for the damping constant «. Clearly this is a sum of 
three terms corresponding to the damping of the free surface (34), the 
side walls (37) and the bottom (40). 

A conventional description of the damping is in terms of the logarithmic 
decrement 6 = 27x/w. ‘The result is 








5 = 6,+6,+5,, (41) 

where 5, = eR (42a) 
v \t/a\(1+(s/k,,, R)? 2k,,<a P 

°s (=) (z){ 1—(s/k,,,R)* sinh 2k,,h f ates, 





v \t/r\ 2k,,R 
5, = (, ) (Z) nh 2k A" (42 c) 
wns sinh 2k,,, A 

‘The case which we have examined in detail is the fundamental mode 
with s= 1, k,,, = 1:841R°'. In figure 2, 6 is plotted as a function of h/R 
for several R; numerical results are given in table 1. 

Three remarks are perhaps in order. 

(a) As would be expected since most of the dissipation takes place at 
the rigid boundaries, the damping is rather insensitive to the shape of the 
cylinder cross-section except in so far as the total area of the walls is affected. 
Thus, the 6 for the circular cylinder is only 18°, less than that for a square 
cylinder of the same cross-sectional area. 

(b) ‘The damping in the body of the liquid (42 a) is small compared to 
the wall damping—particularly for large cylinders. Indeed, in the sense 
of an expansion in powers of //R one should perhaps not retain this term. 
It is readily seen that 5, is O(/?/R*) while 5, and 6, are O(//R). However, 
there are two reasons for keeping it. First, it permits an evaluation of the 
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relative contributions to the damping. ‘Thus, for the smallest of the 
cylinders used, the contribution of this term to the damping was roughly 
10°,. Second, in the limit of an infinitely deep ocean infinite in extent 
this is the only damping term (cf. Lamb (1945), §348). (Of course, then 
k,,, stands for the wave number of whatever wave is being considered and 
is not necessarily zero.) 

(c) The method of calculation adopted is probably not the simplest in 
order to obtain results to the given order in //R. It was used since it is a 
logical first step in a seemingly convergent sequence of successive approxi- 
mations, and permits a clear estimation of the errors remaining after the 
first approximation. 

The comparison of the theoretical prediction (42) with the experiments 
is given below. 


III. ExpeRIMENTAL TESTS 


Measurements were made of the logarithmic decrement as a function 
of h/R for tap water at 16°C in brass cylinders of radius 1-5in. and 3 in., 
and a steel cyliz:der of radius 10in. ‘The brass cylinders were made from 
1/8in. wall extruded brass tubing, the steel cylinder from 3/32in. sheet 
rolled to a 10in. radius and brazed. Bases of 3/8in. thick brass, fitted 
to each of the cylinders by grooving and silver-soldering, were of square 
cross section and cut tangent to the cylinder. 

The amplitude of the surface oscillation was recorded as a function of 
time by means of a sensing element in the liquid, connected through a low- 
frequency amplifier to an Esterline-Angus recorder. The sensing element 
consisted of a 40-gauge wire stretched taut parallel to the axis of the cylinder 
1/16in. from the wall and insulated from the cylinder by a bushing in the 
base. ‘The wire was located accurately on that diameter of the cylinder 
which was perpendicular to one edge of the base. ‘To obtain good linearity 
and stability of response, the wire was connected in series with the cathode 
of the first stage of the amplifier, the conduction path for the vacuum tube 
being completed through the liquid to the grounded cylinder. Since the 
resistance between the wire and the cylinder is proportional (very nearly) 
to the length of wire below the free surface, the change in voltage between 
cathode and grid of the first amplifier stage is proportional to the displacement 
of the free surface of the fluid from equilibrium. ‘To obtain greater sensi- 
tivity for the larger values of 4/R, the wire was extended only 2in. below 
the free surface of the liquid, a rubber band between the wire and the base 
serving to keep the wire taut and to insulate it. ‘Two such wires, each with 
an amplifier and recorder, were used in each tube and were placed accurately 
90° apart. One recorded the motion of the fundamental transverse mode, 
while the other, being placed on the nodal diameter, served as a monitor 
to indicate the purity of the transverse modes. ‘The mechanical arrangement 
is illustrated in figure 1. 

The deflection of the recorder was calibrated in terms of the amplitude 
of oscillation of the fluid by means of a sliding probe fitted with a needle 
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at the lower end and carrying a scale engraved in 0-05in. divisions. ‘The 
probe, located diametrically opposite the sensing element and parallel to, 
and supported from, the wall of the cylinder, was connected electrically in 
series with a generator of 1000 cycle/sec frequency and headphones, the 
final adjustment of the water level being made with an eyedropper. The 
probe was then raised a given distance and the fundamental mode excited 
by rocking the cylinder about the edge of the base. The gain of the 
amplifier was adjusted until the deflection of the recorder equalled a 
convenient number of scale divisions at the time the last intermittent tone 
was heard as the oscillations damped out. By setting the probe at various 
distances above the equilibrium surface, the response of the recording 
system was determined and found to be linear within the accuracy of the 
measurements. 
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Figure 1. Arrangement of the sensing elements and calibration probe. The walls 
of the brass tube are } in. thick, and the base is of 3 in. brass plate. 


Measurements of the logarithmic decrement were made for a range of 
amplitudes, the largest corresponding to a/R approximately equal to 0-2. 
It was difficult to excite larger amplitudes without exciting higher modes. 
However, up to this limit, no deviations from a simple exponential decay 


were observed. 
Figure 2 and table 1 give the results of the measurements for the 1-5 in, 
and 3in. radius tubes. For the number of trials made at each value of 
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h/R the standard deviation of the mean was less than 3% except for 
h/R = 0-25 for the 1-5in. tube; for this it was 9%. The results for the 
10in. radius tube are not included for reasons given below. ‘The build-up 
of a circular mode was essentially eliminated by stressing the cylinder 
slightly to cause the major axis of the ellipse to coincide with the diameter 
containing the sensing element. ‘The eccentricities required to maintain 
the transverse, mode were of the order of a few thousandths of an inch. 
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Figure 2, Comparison of the theoretical and experimental results for the damping 

of surface oscillations in polished right circular cylinders. 


Initial measurements of the decrement with the surfaces unpolished 
gave values of the logarithmic decrement too large by a factor of between 
2 and 3, as shown in figure 3. After polishing the tubes to a mirror finish 
by hand in a direction parallel to the axis of the cylinder the results shown 
in figure 2 were obtained. 
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The bottom surface of the 1-5in. radius cylinder could not be reached 
by hand and was polished in a lathe. This process leaves minute circular 
grooves, which are believed to account for the large values of 5) observed 
for h/R = 0-25 and 0-5, since it is for these values of h/R that the contribution 
of the bottom to the damping is greatest. ‘The decrement for the funda- 
mental mode in this region is actually larger than for the second mode, 
making it more difficult to obtain a pure fundamental. Contamination of 
the inside surface with wax, oil, or even the film left by evaporating alcohol 
increases the decrement by several per cent. 
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Figure 3. Comparison of the theory (solid curve) with measurements in an unpolished 
cylinder. 


The 10in. radius tube, because of the nature of the steel, could not be 
polished to the required degree. The measured decrements were again a 
factor of between 2 and 3 too large, and for this reason the results are not 
tabulated. The shape of the curve obtained was similar to those of figure 2. 
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IV. COMPARISON OF THEORY AND EXPERIMENT 
The overall agreement between theory and experiment (figure 2 and 
table 1) is rather satisfactory and would seem to show that the analysis given 
above is adequate. 














| R(in.) | h/R | f-(sec-) | fo | felfo | 8.x100 | 8 5.- | 8/89 | 5.-/8o 
| 15 | 0-25 2:27 2:26} 1-:005| 7-06 11-8 | 7-08 | 0-598 | 0-600 
me a. 0:50 2:95 3-04 | 0-970 3°45 4-27 | 3-40 | 0-808 | 0-796 
Os, sins 1-0 3:37 3-62 | 0-931 2:47 2-52 | 2-385 | 0-891 | 0-947 
2 1-667 3-46 3-76 | 0-921 2-42 2-28 | 2:32 | 1-06 | 1-02 
Ys 2-0 3-46 3:77 | 0-919} 2-42 2:38 | 2:32 | 1-02 | 0-976 
oy 4-0 3-46 3-80 | 0-911 2-44 2-27 | 2:33 | 1-07 | 1-025 
_ 0-25 1:61 1:73 | 0-932} 4-15 3-89 | 4-00 | 1-065 | 1-025 
i % 0:5 2-08 2:16 | 0-964} 2-02 2:17| 1-98 | 0-933 | 0-914 
i 1:0 2:39 2-47 | 0-968 1-44 1-61 | 1-415 | 0-895 | 0-880 
in 1:5 2-44 2:52 | 0-968 1-40 1-55] 1-38 | 0-904 | 0-891 
= 2-0 2-45 2:55 | 0-962 1-41 1-41 | 1-38 | 1-02 | 0-980 
4-0 2-45 2:56| 0-959} 1-42 1-41 | 1-38 | 1-03 | 0-980 
































Table 1. Comparison of theory and experiment. Surface oscillations in right 
circular cylinders. 


One discrepancy should perhaps be noted. From table 1 one sees that 
the observed frequencies are somewhat higher than those calculated from 
the classical formulae. We believe this can be understood on the basis 
of surface tension effects associated with wetting of the wall. (This does 
not mean the surface tension effects associated with the main part of the 
free surface. ‘This we have seen is much too small.) A crude analysis 
suggests that the correct logarithmic decrement should again be given 
by (42) using as the frequency w the observed frequency. A check of this 
is obtained on noting that one would expect the observed decrement always 
to be greater than the calculated value. (All dirt effects would presumably 
tend to increase the dissipation and hence to increase the observed decre- 
ment.) In table 1 the ratio of theoretical decrement corrected in this way 
(6,.) to the observed (4)) is always (within the standard deviation of the 
measurements) less than, or equal to, unity. 

‘The extreme sensitivity of the experimental results to the condition of 
the walls is quite interesting. Roughnesses whose depth was small com- 
pared to the boundary layer thickness had a remarkably large effect. 


V. CONCLUSION 
It would seem that on the basis of these results one can conclude that an 
essentially correct description of the damping is given by equation (42). 
This is to be understood in the sense of rather ideal situations with very 
smooth walls. Depending on the roughness of the walls, there is an 
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additional factor of between 2 and 4 in the decrement which is approximately 
size independent. ‘This serves to emphasize the extreme care which must 
be used in comparing numerical results of experiments of this kind with 
theoretical predictions. 

Two questions raised by this work are left unanswered. 

(a) Can one understand in detail the effects of very small roughness on 
the damping? 

(6) Can one improve the frequency calculation to obtain closer numerical 
agreement with the observed frequencies ? 
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SUMMARY 

Che finite ditference method due to Lax (1954) is used to solve 
the equations of motion for a cylindrically symmetric flow of a 
compressible fluid. In particular, a converging cylindrical shock 
is found to increase in strength in agreement with the formula of 
Chisnell (1957). ‘The artificial ditfusion introduced by the method 
causes the pressure to remain finite at the axis, but a reflected 
diverging shock is obtained. 


1. INTRODUCTION 

It is well known that a cylindrical shock wave in a compressible fluid 
increases in strength as it converges towards the axis. Some spectacular 
photographs of shocks strengthened in this way have been obtained from 
experiments performed by Kantrowitz (1951). 

‘The problem when the shock has infinite strength has been successfully 
solved by Guderley (1942) and Butler (1954), the solution being singular 
at the axis. For any practical gas there is always some viscosity or heat 
conduction present which will cause the pressure to remain finite 
throughout. 

Neumann & Richtmyer (1950) have described a finite ditference method 
for calculating compressible fluid flows containing shock waves, a feature 
of which is the introduction of an artificial diffusion term. Ina modification 
ot this method by Lax (1954), the coefficient of the diffusion term is 
Av? (2At), where Ax is the width of the space mesh and Af the time interval. 
\lthough this is not the correct viscosity or heat conduction term, the effect 
on the flow will be similar, except in the structure of the diffused shock. 
In these methods the shock is not regarded as an interior boundary, but 
appears as a steep gradient of the variables over a few mesh points. ‘The 
imount of diffusion introduced is very large, as may be seen by considering 
the width of the shock. For air at 300° K, the width is of the order of a 
few mean paths; so that if the flow in one-thousandth of a centimetre of 
air is considered (this distance being covered by 100 mesh points, which 
entails a reasonable length of calculation), the amount of diffusion introduced 
is comparable with what it would be for real viscosity and heat conduction. 

In this paper, calculations by Lax’s method of flows involving converging 
cylindrical shocks are described. It is necessary to give special consideration 
to the flow at the axis, and a reflected diverging shock is obtained in each 
ise. The pressure at the axis remains finite due to the diffusion introduced. 


N 
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‘The converging shock is started by taking initial conditions as in a shock 
tube, namely by considering a cylindrical diaphragm separating two uniform 
regions of a gas at rest with a higher pressure in the outside region, the flow 
being started by the rupture of the diaphragm. Let the ratios of the pressures 
and densities on the two sides of the diaphragm be p* and p*(p* > 1). Ina 
shock tube, if p* = p* (that is, the uniform temperature on the two sides 
is the same), a shock travels into the low-pressure region followed by a 
contact surface, and an expansion fan travels into the high-pressure region. 
By a suitable choice of p* and p*, it is possible to obtain a shock wave and 
expansion fan with no contact surface. 

In this way a converging cylindrical shock of strength 2 is obtained 
(the strength of a shock is defined as the pressure ratio across it). It is 
found that the contact surface does not affect the converging shock, which 
behaves identically for the flow without a contact surface. When stronger 
shocks are considered, the inaccuracy of the method in the vicinity ot the 
contact surface is apparent, and it becomes necessary to eliminate the 
contact surface by having the gas outside the cylindrical diaphragm initially 
at a higher temperature than that inside. For a shock of initial strength $, 
it is found that the inaccuracy of the calculations for the part of the flow 
occupied by the expansion fan also aftects the region of the converging shock 

In each case the calculations were performed with a set of mesh points 
uniformly spaced from the axis to twice the radius of the cylindrical 
diaphragm. ‘The mesh width is taken as 1/64 of the radius. For 140 
time steps of the integration over these 128 points, the time taken on the 
Manchester University Mark I electronic digital computer was 5 hours 
(This included the repetition of each step to provide a check by the 
consistency of the computer.) Each calculation was repeated with a coarser 
mesh of twice this width. One calculation was also performed with half 
this mesh width, taking 16 hours of computing time. ‘The width of the 
shock is proportional to the mesh width. 

‘The strength of the converging shock is found to be in good agreement 
with the formula given by Chisnell (1957). Near the axis the strength 
falls below his predicted value. ‘he amount of entropy left behind the 


outgoing shock is also calculated. 


2. METHOD 
2.1. Differential equations of motion 
The cylindrically symmetric flow of an inviscid compressible gas with 
constant specific heats and without heat conduction satisfies the following 


equations : 


= (rp) + — (rpu) = VU, (1a) 
, ~ j ap 
—(rpu)+ —(rpu*)+r= (), (1 b) 
t cl c) 
— (rE)+ —(rEu+rpu) = 0, (lc) 


} 
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in which r is the distance from the axis, ¢ the time, wu the velocity, p the 
density, p the pressure, and £ is the sum of the internal energy and kinetic 


energy per unit volume: Le. 





where y is the ratio of the specific heats. 
Equations (1) may be made non-dimensional by putting 


r= ror’, Where ro, a constant, is the radius of the diaphragm ; 


p = pop’, Where py, is the initial density of the gas inside the 


diaphragm ; 


~P = pop’, where p, is the initial pressure of the gas inside the 


diaphragm ; 


u=U,u', where uy = (yPo/po)"* 
inside the diaphragm ; 


is the velocity of sound in 


~ 


=f), where? = Fins 


E= pF. 


This gives ( os ey, 
ay (7 p ) (7 pu ) i) 
€ 
— (r'p'u')+ =< (r'p'u?)+y hr’ = = 0 
( ( r 
C ’ va C , v , , , 
<p ( E’)+ = (rE +r'p'u’) = 9, 
( oF 
es h’ ; pamers 
where E nf oe 
ae 


Hereafter these non-dimensional variables will be used and the 


will be omitted. 
It is convenient to introduce new variables a = rp, b = rpu, 
The equations of motion then become 


Ca ch 
: = = 
ct ci 
ch C C 2) 
— + —(bu)+y'?r 2 =, 
ct cr cr 
cc C 
— + —(cu+rpu) =U, 
a + al p 
where a 
| ey 
r 
a 
i=, 
h 
(y- 1) Vy 1) 
p= ——_ - ———__ uy 
r 2 


) 


.2. Finite difference equations 


the gas 


primes 


( rE. 


Consider a set of mesh points r = k Ar, t = t,, where k and n are non- 


negative integers and where Ar is a constant such that A = (Ar) 
integer. Lett,.,—¢, = At,, and let d,, denote d(7, t), where r = r, 


and t= f 
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Let us postpone discussion of the fact that equation (3b) is not in the 
‘conservation form’ desired by Lax (see §2.3), and, otherwise following 
him, replace derivatives by differences as follows. At the point r =r 


t i. replace 





“dy a d, nia (dh, int od, 1) (5 1) 
rt At 
ad d —4¢ 
and = by Pain ia Ln (5b) 
r : 2Ar , # 


We hence obtain explicit equations for the values of the variables a, 6, « 
at time ¢,_, in terms of the variables at time ¢,,, namely 


At, , 
a 1= a Grain) 7A, (1 Pe | eee (6 a) 
iB : 
At. 
Dry st NOx FO ae) 4 Thy 1A un Mtn — ORs Besa 
Yr (Pram —Pkitn)$> (6b) 
At, 
1 Pepin Hepa nd 4 FAp Mo ton p12, Crt esi 


Tay Prawn Mean Ped Pea Mksan) (OC) 
From these equations the values of u, p and p are obtained by equations (4). 
It may be remarked that (i) equations (6) do not apply for k = 0, and 


(ii) the values of variables at points on the two staggered lattices, k +” even 


and k+n odd, are independent of each other. 
2.3. The pressure term 


\s remarked above, the pressure term y !?rcp cr in the momentum 
equation (3 b) is not the derivative of a function of r as is desired for Lax’s 








method. ‘The effect of the pressure is clearer when the ditference equations 





are derived from the conditions of conservation of mass, momentum 
and energy in the space i 
(k— ft) Ar } (k+1)Ar, 
() H i 
where 4 is the azimuthal angle and e <1. In the momentum equation, 


the terms involving the pressure are 


At, ..f A 
Tania Pra Test Praant |. par). 
Possible approximations for [ pdr are 
(Pra + Pra) (7 a) 
2p.» Af, (7b) 
M( Prin t4AP hen + Pra) (7c) 


Phe calculation of a particular flow was done using each of these three 


approximations in turn, and the solutions after 60 time steps were compared. 
It was found that away from the shock the differences between the solutions 
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were negligible, and even at the shock the results for u, p and p differed by 
less than 1°,. In view of this, the approximation (7a), which has the 
advantage that the two staggered lattices are independent, was chosen. 


2.4. The flow at the axts 
At the axis (k = r = 0), we have u = 0 and also a = rp = 0, b = rpu = 0, 
c=rE=0. ‘To obtain the density and pressure it is necessary to resort 
to an approximation different from (6). Several formulae for the density 
and pressure were deduced, and many of these were used in trial calculations. 
The conservation of mass in a cylinder of radius Ar may be written 


al i ee p(r, t) d(r?) + + 2 Arp(Ar, t)u(Ar, t) = 0. 


~ 


ot 


‘Taking the approximation }{p(0,t) + p(Ar, t)}Ar* for this integral, and by 
similarly considering a cylinder of radius 27, we obtain 


At, 1 AI Oo: 
Po,n +1 Pin Tt Ar (43, U3 4 — FP in 41,1): (Sa) 
From consideration of conservation of energy we similarly obtain 
At 
’ MLE 1 lp Z > 
Es 1 E. Ay (ts. Us t 1 Ps. Us , cE, u,, Pi. uy, } (SI ) 
whence Ponsa = (v— 1) Bo.n+1- 


In common with equations (6), equations (8) use a staggered mesh. 

When trial calculations were made, it was found that in the pressure 
and density at the axis, oscillations developed and rapidly increased in 
amplitude. (The oscillations for the two independent staggered lattices 
usually differed in phase and amplitude.) It was found that this numerical 
instability appeared to be eliminated by using the approximation (7 c) at 
the pointk = 1. Thus, for k = 1, equation (6b) is replaced by 


Dy ns = $00n + a | — by, Us, +y**Ar(4 Po, 4 +p, 3 Pon)}- (9) 


[n all calculations done with this modification to equations (6), oscillations 
at the axis did not appear. However, the term in p, ,, makes the two staggered 
lattices no longer independent of each other. 


2.5. Initial conditions 

‘To calculate the flow near the axis it was found necessary to perform 
the integration simultaneously for the two staggered lattices. (That is, the 
calculation is performed at all the points k = 0, 1,2... at each step in time 
n=,1,2....) This adds importance to the precise initial values of the 
variables in the neighbourhood of the diaphragm. At the diaphragm k = K, 
the variables are chosen so that ag 9, Dg», Cx are the respective averages for 
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k= K+1 of a), b).9, ¢;9. Thus the initial conditions are taken as 
uj.9 = 0, for all k, 
Pro=1,pPpo=1, for k<K, 


Pk.o = pr, Pro = a, for k K, 


vhere p* and p* are constants. 


2.6. The time interval 


In his paper Lax mentioned the Courant-Friedrichs-Lewy criterion as 
i necessary condition for the stability of this numerical method. In this 
problem the largest velocity is that of the shock, so that this necessary 
condition for stability is 

At, /A l (velocity of shock wave). 
Also, the width of the shock depends on the value of A¢,,/Ar, and is smallest 
vhen this is as large as possible. 

\ series of trial calculations were made choosing At. so that 

\t,, Ar = A/(velocity of shock wave), 

where A assumed various constant values. For a weak shock the method 
Mf calculation was unstable with dA = 0-9, but with 4 < 0-85 it appeared 
to be stable. (Any instability first showed at the rear of the shock.) Hence 
the time interval was chosen as above with A = 0-85. For stronger shock 
waves this choice of time interval also gave instability, and for the strongest 
shock (initial strength 8) it was found necessary to further reduce the time 
interval by taking 4 = 0-75. 

Throughout, the time interval was chosen for the next step from the 
shock velocity in the previous step, which does not involve any iteration. 
Uhis shock velocity was calculated from the distance moved by the shock, 
the position of the shock () R) being determined from the pressure 
distribution. ‘This was done by linearly interpolating in the table of p, 
against 7, to find the point where p is the average of the pressures behind 
and in front of the shock. The pressure behind the shock is taken as the 
local maximum value of p;.,,.. For the shock converging towards the axis, 
the pressure in front of the shock was taken as p,,,, the pressure at the axis. 
For the diverging shock, the pressure in front of the shock was defined as 
the pressure at the point where p(y — Av, t,,)—p(,t,,) = B, where suitable 
values for the constant B were found to be: 0-1 for a weak shock of strength 2 
initially; 0-2 for a shock of initial strength 4+; 0-4 for a strong shock of 
strength 8 initially. 

When near the axis, the position of the shock is not identified sufficiently 
.ccurately. ‘The procedure adopted was to keep the time interval constant 
at the value it has obtained as the axis is approached. When At, was kept 
it this small value, the reflected shock was so diffused that it was impossible 
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to recognize it. ‘lhe outgoing shock travels at a lower velocity, thus 
permitting a larger time interval. ‘Therefore, after a few steps, a larger 
constant time interval At, = C was taken. The value of C was chosen so 
that when it was possible to adjust the time interval from the velocity of 
the outgoing shock, the value of At,, was then slightly greater than C. 

This method of adjusting At, was found to be sufficient to prevent the 
appearance of numerical instability due to the use of too large a time interval, 
and to avoid excessive spreading of the shock due to taking At, much smaller 
than the permissible maximum. 


3.1. Initial pressure and density ratio 4 

For a shock tube, the initial conditions specified by p* = p* = 4 give 
a shock of strength 1-93, a contact surface and an expansion fan. The 
solution is obtained for a cylindrical diaphragm with these initial conditions, 
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Figure 1. Pressure vs radius at intervals 0-2 of time for a flow initiated by a cylindrical 
diaphragm with initial pressure and density ratios 4. The converging shock 
of initial strength 1-93 increases in strength and reaches the axis at ¢ = 0°66, 
and a reflected shock is obtained. An expansion fan travels out from the 
diaphragm. The width of the space mesh used is 0-008 units where the 
radius of the diaphragm is 1 unit, the pressure ahead of the converging shock 
being the unit of pressure. 
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a converging cylindrical shock initially of strength 1-93 being obtained. 
The mesh width Ar= 1/128 is used. 

In figure 1 the pressure distribution is shown at intervals 0-2 of time 
(each representing about 40 steps of the integration). Apart from the 
first few steps of the integration, the shock appears as a rapid variation 
in p over about six points of the space mesh. As f increases and the shock 








Figure 2. Velocity vs radius at intervals 0-2 of time for a flow initiated by a cylindrical 
diaphragm with initial pressure and density ratios 4. The unit of velocity 
is the velocity of sound in the undisturbed gas ahead of the converging shock. 
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Figure 3. Density vs radius at intervals 0-2 of time for a flow initiated by a cylindrical 


diaphragm with initial pressure and density ratios 4. The contact surface is 
spread over a large number of mesh points. The density is unity ahead 
of the converging shock. 
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moves in towards the axis, the strength of the shock is seen to increase. 
After the shock has passed, the pressure at any fixed point behind it continues 
to rise. When the shock reaches the axis, the pressure there rises to a high 
but finite value (figure 5), and a reflected diverging shock appears. The 
pressure at a fixed point behind the reflected shock decreases as ¢ increases. 
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Fisure 4. Temperature vs radius at intervals 0-2 of time for a flow initiated by a 
cylindrical diaphragm with initial pressure and density ratios 4. The diverging 
shock leaves a region of heated gas between the axis and the contact surface. 
The initial temperature of the gas is unity. 


In figure 2 the velocity of the gas in this flow is seen to behave similarly, 
except of course that the converging shock decreases the gas velocity from 
zero to a negative (inwards) value. The diverging shock increases the 
gas velocity to a small positive (outwards) value. Behind the converging 
shock, the velocity at a fixed point also increases with time and decreases 
when the diverging shock has passed. 

The graphs of the density (figure 3) are similar to those of the pressure, 
but the rise in density across the shock is smaller, corresponding to a 
temperature increase. In the distributions of both density and temperature 
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(figure +), the contact surface appears. Unlike the shock, this is far trom 
being a discontinuity in these variables; it is a gradual change spread over 
in increasing number of mesh points. The contact surface moves inwards 
behind the converging shock, and is traversed by the diverging shock. 
\t t = 1-4, the shock has more or less passed the smudged contact surface, 
leaving a region of high temperature between the axis and contact surface 
(see 3.5) 
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Figure 5. Pressure, density and velocity vs time at the fixed point r — 0:375. At the 
axis r — Q, the two rapid changes in p and p when the converging and diverging 
shocks pass are merged. ‘The velocity is zero at the axis. 





The variation with time of the variables at fixed points r = Q and 
r = 0-375 is shown in figure 5. Although the general features resemble 
the solution given by Guderley (1942), the magnitudes of the variables 
are considerably smaller since the condition of infinite shock strength is 
not fulfilled. 
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3.2. Stronger shock: heated external region 

With the above initial conditions (p* = p*), a contact surface is produced 
which cannot be satisfactorily treated by this method of Lax (1954). With 
the shock wave, where the density, velocity and pressure decrease in the 
direction of motion, there is the tendency of the flow to make the wave 
\arrower, Which balances the spreading effect of ditfusion. With a contact 
surface, there is no such balancing effect and diffusion continues unimpeded. 
When a higher pressure ratio was used, the spreading of the contact surface 
increased and appreciably affected the shock. In order to eliminate the 
indesirable contact surface, values of p* and p* were chosen such that the 
shock-tube equations are satisfied with a contact surface of zero strength. 
For example, the values p* = 3-52, p* = 2-44 satisfy the shock-tube 
equation with no contact surface and a shock of strength 1-93 (which 1s 
the same strength as that given by p* = p* = 4). With these initial 
-onditions, the converging cylindrical shock behaves identically. Differences 
in the diverging shock only occur when in the former case (p* = p* = 4) 
the contact surface is traversed. On passing into a colder region, the shock 
travels more slowly. Similarly, the outgoing expansion fan travels more 
juickly in the latter case when the external region is at a higher 
temperature. 

By initially heating the external region to eliminate the contact surface, 
it is possible to obtain a stronger converging shock, and results were obtained 
tor flows with shocks of initial strengths 4 and 8. A mesh width Ar = 1/64 
vas used. ‘The distributions of pressure, velocity, density and temperature 
t intervals 0-1 of time are shown for the latter example in figures 6 to 9. 
The general features of the behaviour of the variables are the same as for 
the weak shock (strength 1-93 initially). As predicted by Chisnell (see 
$3.3), the stronger shock increases in strength more rapidly than the weak 
shock as it approaches the axis. 

From the pressure distributions (or perhaps more clearly from the 
velocity distributions), it will be seen that for the stronger shock flow the 
expansion fan has increased diffusion. (‘This is partly but not entirely due 

the use of a coarser space mesh.) ‘This diffusion is particularly notice- 
le in the first few steps of the integration. It is not unreasonable to 
suspect that the diffusion from the expansion fan may affect the converging 


s! ock. 


3.3. Strength of converging shock 
‘he position and strength of the shock at any instant of time are deduced 

om the pressure vs radius distributions such as are shown in figures 1 and 6 
(see $2.6). The strength of the shock is taken as the maximum value of the 
pressure in the p/r distributions. Because of the spreading of the shock 
over several mesh points, the strength of the shock will be underestimated 
when the axis is approached, and the position of the shock is also possibly 
in error by the order of one mesh length. 
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Figure 6. Pressure radius at intervals 0-1 of time for a flow initiated by a cylindrical 
diaphragm with initial pressure ratio 38:2 and density ratio 10°7. The 
converging shock begins at strength 8 and reaches the axis at f 0-31 
\luch greater rises in pressure are obtained with this stronger shock The 
width of the space mesh is 0-016 and the pressure ahead of the shock is unity. 








Figure 7. Velocitv vs radius at intervals 0-1 of time for a fow with a converging 
evlindrical shock of initial strength 8. The velocity of sound ahead the 
converging shock is unity. 
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Figure 8. Density vs radius at intervals 0-1 of time for a flow with a converging 
cylindrical shock of initial strength 8. The density ahead of the converging 
shock is unity. 
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Figure 9. Temperature vs radius at intervals 0-1 of time for a flow with a converging 


cylindrical shock of initial strength 8. Initially the temperature between the 
axis and diaphragm is unity, and beyond the diaphragm 1s 3°57. 
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In view of the crudeness of this part of the calculation, it is perha 


remarkable that the results obt 


tf 
ained should fit a smooth curve of z cs R, 
where z is the shock strength and R js the radius of the shock (figures 10, 1 
The agreement with the formula (obtained by Chisnell (1957)) for t 
strength of a converging cylindrical shock is extremely good. ‘T} 
divergence from | 


1) 
i 


us predicted strength only appears near the axis and 
decreases with the mesh width. 
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Figure 10. Shock strength radius for a converging cylindrical shock of initial 


Strength 1-93. The agreement with Cl 


usnell’s formula is better when a finer 
mesh 1s used 


lor a stronger shock the agreement is not so outstandingly good. The 
strength obtained is greater than Chisnell’s. This deviation is attributed 
to the effect of diffusion from the expansion fan (see § 3.2 
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Figure 11. Shock strength vs radius for a converging cylindrical shock of initial 
strength 8. Errors arise from the diffusion of the expansion fan. 


3.4. Variation of gas constant y 

The results described above are for a gas with y= 7/5. The calculations 
were repeated for the flow with initial conditions specified by p* = p* = 4, 
with y = 5/3. The pressure distributions were very similar to the former 
flow, the initial shock being just slightly weaker. The variation of shock 
strength with distance from the axis was again in good agreement witl 
Chisnell’s theory. A larger density change across the contact surface was 
obtained. The temperature near the axis when the diverging shock had 
traversed the contact surface was much higher, being 1-95. The larger 
temperature variation was accompanied by a smaller gas velocity. 


3.5. Entropy increase due to cylindrical shock 


When the reflected diverging shock has travelled out, a region of gas at 
high temperature and almost at rest is left between the axis and the contact 
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surface (figures 2 and 4). ‘The |.» temperature indicates that a con- 
siderable amount of energy has be. »: «ssi pated by the shock. For a pertect 
gas without any diffusion, the flow ntropic except on passing through 
the shock. ‘There is no change in the entropy per unit mass in flow through 
the expansion fan, and there is no fiow through the contact surface, at which, 
however, there is a discontinuity in entropy. Hence, when the diverging 
shock has passed, the total entropy of the gas between the axis and the 
contact surface will be constant (assuming no subsequent shocks pass 
through this region). ‘This heated gas is the gas which was initially between 
the axis and the diaphragm, the shock having passed through it on both its 
inward and outward journeys. 

‘The entropy increase produced by the shock is calculated as follows. 
lhe contact surface r = x is located by the fact that the mass of gas between 
the axis and contact surface is constant, that is, by finding x such that 

J 
2nprdr = | 2ardr =z. 
0 


‘The total entropy per unit length in this gas is 


2zc,{log p— y log p)pr dr, 
o 
where c, is the specific heat at constant volume. ‘The integrals were evaluated 
by Simpson’s rule, and, for the flow with initial conditions specified by 
p* = p* = 4, gave large negative values for this entropy which also decrease 


with time. ‘The initial entropy of this gas is zero, and the entropy per unit 
mass behind the initial shock is 0-O11c, and is greater tor stronger shocks. 
However, the entropy per unit mass in the gas initially outside the diaphragm 
is —0-555c,, and diffusion of entropy across the contact surface exceeds 
the production of entropy by the shock. 

For the How with the same initial shock strength but with the gas 
initially heated to give a contact surface of zero strength, the entropy per 
unit mass outside the diaphragm is 0-011¢, which is the same as the entropy 
behind the initial shock. Diffusion across the contact surface will be 
reduced, making this a more suitable flow from which to calculate the 
entropy produced by the shock. ‘The total entropy per unit length left 
between the axis and contact surface is found to be 0-077c,. This may be 
compared with a plane shock of strength 1-93 reflected normally by a plane 
solid wall. ‘The increase in the total entropy in a mass 7 of gas is 0-063. 


for a plane shock reflection. 
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Upward ‘ falling’ jets and surface tension 


By JOSEPH B. KELLER and MORTIMER L. WEITZ 
Institute of Mathematical Sciences, New York University 


(Received 1 November 1956) 


According to the simple hydraulic theory of jets, each particle of a jet 
moves independently along a parabolic trajectory. Therefore a steady jet 
has a parabolic shape. We wish to consider how these results are modified 
by surface tension. For simplicity we will consider a two-dimensional jet 
of incompressible fluid. 

The hydraulic theory is based upon the two assumptions, that the 
velocity is constant on each cross-section, and that the pressure is constant 
throughout the jet. ‘Then the jet can be completely described in terms of 
its centre line (x, ¢), its vertical width A(x, t), its horizontal velocity u(x, t) 
and its vertical velocity z(x,¢). These four functions satisfy the following 
four differential equations, the first of which is a conservation of mass 
equation, the second a kinematic relation, and the last two are momentum 
equations. 


ph, + p(uh),. = 0, (1) 

yy, tuy, =v, (2) 
ph(u,+ uu) 2Ty, ¥24A1+y%)*?, (3) 
ph(v,+uv,) = 2Ty,,(1 + y2)-3? — pgh. (4) 


In these equations p denotes the density of the liquid, 7’ denotes the tension 
of the surface, and g denotes the acceleration of gravity. 

Let us seek a steady (1.e. time-independent) solution of these equations. 
From (1) we obtain uh = m, where m is the constant flux of mass through 
any cross-section of the jet. Equation (2) yields v = uy,, while (3) and (4) 


r 





become 
—2T E 
u, = oF +h (5) 
m : 
2F Pee 
(uy,). = om deel + y7) 3? — pur", (6) 


Upon integrating (5) and denoting u(0) by uy, we obtain 
_ 

u(x) = y+ SCL +98) 12 (1+ y90))79) (7) 

Now on using (7) to eliminate u(x) from (6), we have 
[l+a(l+y?)2]y,.=—B. . (5) 
Here the constants x and § are defined by 
ar 
muy — 2T(1 + y2(0))-1?’ 
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om 
6 —_—_— =) ( ] () ) 


[mu —27(1 + y%(0)) 12 





Integration of (8) yields 


y,tasinh ly, —Px+y(0)+asinh“! y,(0). (11) 

The last equation is a first order equation for y(x) involving three 
parameters. ‘To reduce the number of parameters we set y,(0) = 0, since 
the resulting curve will apply to any initial slope if read from that slope on. 


Q Q 


Next, we introduce the new variables Bx = € and By = 7. In terms of 


these variables (11) becomes 


iget asinh™! yg = —€. (12) 


This equation involves only the one parameter x = 27)/(muy—2T) and is 
therefore convenient for numerical integration. 

When «= 0, (12) yields the expected parabola 7 sé*. In the 
figure graphs of the solution of the above equation are shown for « = 0, 1 
and —2. When the €, 7 variables are used, the curves with « > 0 lie above 
the parabola obtained for z = 0. On the other hand, when the x, y variables 
are used these curves lie below the parabola. ‘Thus, as one expects, in the 
physical plane the jets with surface tension lie below the parabola, provided 
x4 > 0, or mu,/2T > 1. 

However, when x — 0, or mu,/27 < 1, the jet rises instead of falling, 
even though it is initially projected horizontally (see the figure). It continues 
rising until all the initial kinetic energy is converted into potential and 
surface energy. ‘hen the curvature and thickness become infinite and 
the theory fails. Presumably it would spill down before this point is 
reached. ‘This phenomenon of a rising jet can occur only in a slow thin 
jet, since the above inequality requires that the kinetic energy }phu? must 
be less than the surface energy 7. In very crude experiments we have been 
unable to observe this behaviour. We believe this is because the rising 
flow is unstable at the necessarily low speed. Instead, a ‘teapot-like’ 
flow, in which the liquid runs along the lower surface of the spout, seems to 
be stable. 

In order to understand the surprising phenomenon of a jet ‘falling’ 
upward, let us consider the simple problem of a falling body of mass M. 
If \(t) denotes the height of the body measured positive upwards, then the 
equation of motion is My, Mg. From this equation we see that v, 
is negative. If we want the body to fall faster, we may push it down with 
a force k*y, proportional to the acceleration. ‘Then the equation of motion 
becomes 

Mv, = —Me+k*y,, 
and we now have Vy = —Me/(M —R?*). 


Thus as Rk? increases from zero, the acceleration becomes more negative, 
until k? V/, when there is no solution. For k? > M, yy is positive so the 
body falls upward. We may describe this example by saying that negative 
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inertial mass of amount — k” has been added to the body. As more negative 
inertial mass is added, the body’s downward acceleration continually 
increases until the total inertial mass M—k? becomes negative. ‘Then the 


body falls upward. 











€ 


Steady jets falling or rising under the influences of gravity and surface tension 


The parabola « = 0 occurs when surface tension is absent. For moderate 
surface tension (« ~. 0) a curve such as that shown for a 1 occurs. In 
terms of the x,y variables this curve lies below the parabola. However, 
for larger surface tension (a — 0) a curve like that shown for « 2 results. 
This curve represents an upward falling jet. The curves for a 1 and 
¥ 2 were obtained by numerical integration of equation (12). 


This example is pertinent to the jet problem. ‘lo see this, note that 
near the nozzle the x-coordinate of a particle is approximately u,t. ‘The 
surface tension force is approximately 27y,. = 2Ty,/u>. Thus the above 
discussion applies to a particle of the jet with k? = 27/u5. We see that 
surface tension has the effect of adding a negative inertial mass to each 
particle without changing its gravitational mass. 

These results appeared in ‘‘ Thin unsteady heavy jets”’, by J. B. Keller 
and M. L. Weitz, Report IMM-NYU 186, Institute of Mathematical 
Sciences, New York University, 19 December 1952. ‘They were also 
presented at the Ninth International Congress for Applied Mechanics, 
Brussels, September 1956. 
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REVIEWS 


Momentum Transfer in Fluids, by W. H. Corcoran, J. B. OPrELL 
and B. H. Sacer. New York: Academic Press Inc., 1956. 394 pp. 


72s. or $9.00. 


‘Towards the end of his monumental work on Hydrodynamics, Lamb 
says, ‘‘ it remains only to call attention to the chief outstanding difficulty 
of our subject ’’, and goes on to explain that he is referring to turbulence. 
‘That sentence appeared first in the 1895 edition, and it has remained in all 
subsequent editions. Having called attention to the * difficulty ”’, Lamb 
let the matter rest there and said very little more. Lamb was fortunate, of 
course, in that he had a good story to tell about numerous other parts of 
the subject, and nobody minds the omission. Other, later, writers of 
textbooks on fluid mechanics for the student or the practicing engineer are 
not so well placed, and are obliged to do more than “ call attention to the 
difficulty”. Some help must be provided to the aeronautical engineer 
who wishes to calculate skin friction, to the meteorologist who wishes to 
predict atmospheric diffusion rates, and to the chemicai engineer who 
wishes to relate mass transfer in a model experiment to that in a full-scale 
plant. Modern technology needs help in describing and analysing turb- 
ulent flow, and cannot wait for scientists to understand its mysteries. Nor 
is it sufficient to provide a collection of data in the style of a handbook, 
since the variety of types and effects of turbulent flows is too large to allow 
cataloguing. Some kind of rational analysis of turbulence is demanded, 
in such a form that the student or engineer can apply it with understanding 
and confidence in many different situations. 

Kaced with this need, the modern writer of a textbook on turbulent 
Hows of one kind or another can proceed in either of two different ways. 
The first way, which has been commonly adopted, is to present the older 
analysis usually described as ‘ mixing-length theory’. ‘This analysis 
serves quite well as a framework for an empirical description of turbulent 
motion, but as a theory it has lost vitality and has been at a standstill for 
20 years. It seemed at one stage to have considerable promise, since it was 
able to provide reasonably accurate predictions of mean velocity profiles 
in wakes, jets, mixing regions, and other cases of nearly unidirectional 
mean flows, with only a small number of assumptions; however, later 
work suggested that this was due in large measure to the insensitivity of 
the profiles to the assumptions used, and that, as a consequence of the 
strong geometrical restrictions on the possible form of the profile in such 
cases, almost any theory can achieve some success. ‘The engineer is 
especially interested in transfer rates, and here ‘ mixing-length theory ’ is 
much less successful in its predictions; in this sphere it has value only as 
a language for the description of what has already been measured. 

It is all to easy to adopt a superior attitude to the old ideas of transfer 


of momentum, vorticity and heat by a mixing process. They have been 
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fair game for attack by many of us for a number of years, but no one has 
had much success in producing an alternative line. In 1950 Schlichting 
wrote, in the preface of his book Grenzschlicht-Theorie, ‘‘ It is true that 
according to present views these theories possess a number of shortcomings 
but nothing superior has so far been devised to take their place, nothing that 
is, which is useful to the engineer. ... statistical theories of turbulence have 
not attained any practical significance for engineers”. I do not think 
the second part of this quotation was valid in 1950, and I am sure it is not 
valid now; however, there will be general agreement with the first part, 
and therein lies the rub. ‘The modern or so-called statistical theories of 
turbulence (not entirely an appropriate description since there is no parti- 
cular theory or set of assumptions or hypotheses underlying the modern 
work; it is characterized by nothing more than a recognition of the random 
character of turbulent velocity fluctuations and of the necessity to describe 
turbulence statistically) have not yet fulfilled the promise they were 
believed to have. ‘The first aims of the ‘statistical theory’, to develop a 
mathematical framework for the description of a random three-dimensional 
continuous velocity field, and to use the equations of motion and of cont- 
inuity to obtain relations between different elements of this description, 
have been achieved, and many people believe that this prepares the way for 
the development of theories that will yield useful results for the engineer. 
So far, that has not happened to any appreciable extent, although recent 
advances in understanding of the mechanism of turbulence, as described 
by ‘Townsend in his book The Structure of Turbulent Shear Flow, give 
renewed hope that the engineer will not have to wait very much longer. 
statistical theory ’, whatever its 


There is of course no going back on the 
lack of success in hard practical terms, since in essence it is only a descriptive 
framework, whose use is just as inevitable in turbulence as it is in any other 
subject involving random quantities, such as molecular motion in gases. 
‘The second possible course open to the writer of a textbook about 
turbulent flow is to grasp the nettle firmly, to do the best he can with 
modern ‘ statistical theory ’ and to refuse the doubtful help of discredited 
or unreliable theories. ‘This course has not yet been pursued whole- 
hearted in a textbook, and we cannot readily judge whether it would be 
successful in meeting the needs of engineers. It may not be quite as 
difficult as at first it seems. By making the widest possible use of dimen- 
sional analysis (which reproduces many of the intermediate results of the 
mixing-length theories), and by introducing inferences from previous 
measurements in very general form, particularly those concerning firstly the 
region of nearly uniform stress near a rigid wall and secondly the region 
of a turbulent flow that is not affected directly by the presence of a rigid 
wall (the ‘law of the wall’ and the ‘ law of the wake ’, discussed recently 
by Coles in this journal, are examples of such empirical generalizations), 
all this being done within the framework provided by the ‘statistical 
theory’ and interpreted in the light of the now considerable body of 
knowledge about the mechanics of turbulence, it might be possible to 
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produce a satisfactory textbook book about turbulence. It is certainly 
worth trying, in view of the potential value of such a book. 

The book under notice is by three chemical engineers at California 
Institute of ‘technology and, according to the preface, is intended to 
describe a part of fluid mechanics that is ‘‘ of special interest to chemical 
engineers’. ‘They take this to include the fundamentals of viscous fluid 
flow, boundary layers, and turbulence, with emphasis on cases of flow 
through pipes and channels. ‘That three teachers of chemical engineering 
should feel impelled to write a textbook, mostly about turbulent flow, 
for use by advanced students and practicing engineers, suggests that the 
writings of turbulence specialists have not been found very helpful. Such 
a venture would seem to demand courage, but in fact a reading of the book 
suggests that it is the courage of the blind who do not know what pitfalls 
lie about them. At a brisk pace they cover a wide range of topics of 
alarmingly different degrees of difficulty, and spend little or no time on 
discussion of the pros and cons of different ideas. When an hypothesis is 
introduced, the responsibility is pushed back on its originator; von Karman 
said this, Prandtl said that, this is Schlichting’s result, that is Pohlhausen’s, 
and soon. ‘This is the handbook manner, which is very useful in its place, 
but which leads here to a disturbing lack of discrimination in the selection 
and presentation of the available information about turbulence. 

‘The authors seem to have been uneasily aware of the above two possible 
points of view about a textbook treatment of turbulence, and their answer 
is to use both. Not in the same chapter, for that is scarcely possible, but 
in different parts of the book. Of the five large sections on turbulent flow 
(four chapters and an appendix), three are substantially about the pre-war 
theories and types of measurement, and say little that is not in Modern 
Developments in Fluid Dynamics (1938) or Bakhmetett’s book, The Mechanics 
of Turbulent Flow (1941); the other two use the language and notions of the 
‘modern approach’ to turbulent flow. ‘These two divisions of the book are 
not integrated, or even related, and one suspects that the authors were playing 
for safety; if the reader wants the simple relations and analysis provided by 
transport theories, the book provides that, and, if tomorrow’s research 
should suddenly make statistical theories useful, all the reader has to do 
is to jump a chapter or two. Of the two divisions, that concerned with the 
‘modern approach ’ is the less satisfactory. ‘The authors evidently have 
great respect for the complicated game played with tensors by theoreticians, 
and their deference has led them to include a large appendix containing 
long derivations of formulae (for example, for triple correlations 1n isotropic 
turbulence) of which they make no use and about which they make no 
comment. ‘This is no more than a genuflexion before the altar of ‘ high- 
brow ’ analysis, and the book would be both more honest and less confusing 
without it. The chapter on the ‘modern approach’ to turbulence is 
similarly characterized by enthusiasm rather than care, and by an absence 


of any comment which would help a reader to judge the value of the 
ditferent concepts and ideas. Even the experimental data reproduced 
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here is of doubtful value, since the quantities concerned, such as the 
kinetic energy spectrum, are not explained adequately. No; so far as 
current ideas about turbulence are concerned, this book is not likely to 
help the student or engineer; it is unlikely either to make available 
to him the small stock of practically useful results or to give him an 
understanding of what the current ideas are about. 

Such contribution as the book does make, and it cannot be reckoned 
as large, is to gather together the data about flow in pipes, channels and 
boundary layers, including the effects of wall roughness and wall heating, 
and to present this in the way that has been familiar to aeronautical and 
mechanical engineers during the past 25 years. A book with this aim might 
well be very useful, since most of the existing accounts, such as Modern 
Developments in Fluid Dynamics, perhaps give more discussion than is 
useful to a person who wishes only to know what data is available and when 
it is applicable. I think that this was what the authors set out to do, 
but I cannot feel that they have succeeded. ‘They seem to have paid little 
attention to the p/an of their work or to the need for systematic presentation 
of material. Chapters 2, 3 and 6 are entitled ‘* Some simple properties 
of turbulent flow”, ‘“‘Some macroscopic characteristics of turbulent 
flow’, and “‘ Some properties of turbulence ”’ ; how can anyone find his 
way about a book that provides such vague directions ? Within each 
chapter the arrangement of sections has little perceptible order or of the 
gradual development of a theme. ‘The vorticity transfer hypothesis is 
introduced and, in order to put the reader in the picture, the definition of 
vorticity has to be provided in a footnote; the no-slip condition at a rigid 
wall is given in the section entitled ‘“‘ Boundary conditions for the equation 
of continuity ’’; a section entitled “‘ Measurement of the physical nature 
of turbulence ”’ (whatever that may mean) begins with talk about velocity 
fluctuations and then for no apparent reason proceeds to devote most of its 
length to a discussion of the width of the mean temperature profile behind 
a hot body in a turbulent stream; and many other examples of an almost 
random arrangement of material could be given. ‘The authors had the 
advantage of writing in 1956, with many more published papers to call 
on for help, but even so I doubt if their book is clearer and more inform- 
ative than the many older accounts of the subject. 

The book is also marred by slipshod writing. Non sequiturs, clumsy 
constructions, ill-chosen words, and even grammatical errors, are numerous 
in all parts of the book. We are accustomed nowadays to an unadorned, 
rather blank, arid, style of writing in scientific books, and for some reason 
books for engineers seem to have gone farthest in this direction, but surely 
readers may still expect a certain minimum standard of care and respect for 
the rules of composition. ‘To reach this standard dogs not require literary 
skill and practice; it requires only that one should be willing to take time 
and trouble, and in particular to go over one’s work several times. Is it 
conceivable that each of the three authors read and approved of the 
sentence (p. 114) ‘‘ The restriction homogeneous is used so that the equations 
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to be developed are mathematically correct inasmuch as the analyses 
will be confined to regions where all functions of interest are continuous 
and differentiable” ? One sees what the authors mean in cases where 
precision of thought is not needed, but there are places where their sloppy 
writing definitely impedes understanding. I can only guess the meanings 
of ‘‘ The intensive properties of the fluid will be considered to be inde- 
pendent of state’’ (p. 11), of the symbol j defined as “ Friction per unit 
time in differential length of total conduit, ft. lb./sec.” (p. 44), and of 
‘Steady laminar boundary layers are actually of infinite thickness but it 
is convenient to consider finite thicknesses of the boundary layer for 
purposes of practical calculations’ (p. 240). The effect on a measuring 
instrument of its velocity relative to the fluid is described as if it were a 
property of the coordinate system (p. 7), boundary layer separation is 
entangled with stability (p. 265), and § VII-9 is a thorough muddle of 
the two cases of a time-dependent laminar boundary layer on an infinite flat 
plate and of a steady boundary layer on a semi-infinite flat plate. Many 
other examples could be given, and remarks in similar vein could be made 
about the clumsiness of the more mathematical arguments. ‘This is 
‘ first-draft ’ writing and it is sheer irresponsibility to put it into print. 

It is of course always possible to learn something trom a new book. | 
myself picked up one little wrinkle, which I have stored away mentally 
for use on occasions when I wish to impress students with my progressive 
outlook; it is that when one is establishing the equation of continuity, it 
is nowadays desirable to exclude the possibility that nuclear reactions are 
going on in the volume element under consideration. 


G. K. BATCHELOR 




















JOURNAL OF 
FLUID MECHANICS 


Orders originating 


in 
UNITED STATES OF AMERICA 


and 


CANADA 


should be sent to the sole agents 


ACADEMIC PRess INC, 
111 Fifth Avenue, 
New York, 3, 
N.Y., U.S.A. 














JOURNAL OF FLUID MECHANICS 
Volume 2, Part 2 March 1957 


CONTENTS 


An experimental investigation of heat transfer effects on 
boundary layer separation in supersonic flow 
by G. E. GADD 


The Saha equation and the adiabatic exponent in shock 
wave calculations 
by RALPH A. ALPHER 


Buoyant plumes in a moist atmosphere 
by B. R. MORTON 


A new approach to problems of shock dynamics. Part I 
Two-dimensional problems 
by G. B. WHITHAM 


Damping of surface waves in an incompressible liquid 
by K. M. CASE and W. C. PARKINSON 


A numerical method for a converging cylindrical shock 
by R. B. PAYNE 


Upward ‘ falling ’ jets and surface tension 
by JOSEPH B. KELLER and MORTIMER L. WEITZ 


Reviews .. 











